OSM report 2022-2023
Marissa Dyck
2024-05-08
The first two chunks of this r markdown file after the r setup allow for plot zooming, but it also means that the html file must be opened in a browser to view the document properly. When it knits in RStudio the preview will appear empty but the html when opened in a browser will have all the info and you can click on each plot to Zoom in on it.
Before you begin
Notes
A few notes about this script.
If you are running this with the 2022-2023 data make sure you download the whole (OSM_2022-2023 GitHub repository)[https://github.com/ACMElabUvic/OSM_2022-2023] from the ACMElabUvic GitHub. This will ensure you have all the files, data, and proper folder structure you will need to run this code and associated analyses.
Also make sure you open RStudio through the R project (OSM_2022-2023.Rproj) this will automatically set your working directory to the correct place (wherever you saved the repository) and ensure you don’t have to change the file paths for some of the data.
Lastly, if you are looking to adapt this code for a future year of data, you will want to ensure you have run the 2_ACME_landscape_covariate_exploration_script.Rmd with your data as there is much data formatting, cleaning, and restructuring that has to be done before this code will work. Helpful note: The files are numbered in the order they are used for this analysis.
If you have question please email the most recent author, currently
Marissa A. Dyck
Postdoctoral research fellow
University of Victoria
School of Environmental Studies
Email: marissadyck17@gmail.com
(update/add authors as needed)
R and RStudio
Before starting you should ensure you have the latest version of R and RStudio downloaded. This code was generated under R version 4.2.3 and with RStudio version 2024.04.2+764.
You can download R and RStudio HERE
R markdown
This script is written in R markdown and thus uses a mix of coding markup languages and R. If you are planning to run this script with new data or make any modifications you will want to be familiar with some basics of R markdown.
Below is an R markdown cheatsheet to help you get started,
R
markdown cheatsheet
Install packages
If you don’t already have the following packages installed, use the code below to install them. *NOTE this will not run automatically as eval=FALSE is included in the chunk setup (i.e. I don’t want it to run every time I run this code since I have the packages installed).
install.packages('tidyverse')
install.packages('ggpubr')
install.packages('corrplot')
install.packages('Hmisc')
install.packages('glmmTMB')
install.packages('MuMIn')
install.packages('TMB', type = 'source')
install.packages('rphylopic')
install.packages('broom')Load libraries
Then load the packages to your library so they are usable for this session.
library(tidyverse) # data tidying, visualization, and much more; this will load all tidyverse packages, can see complete list using tidyverse_packages()
library(ggpubr) # make modifications to plot for publication (arrange plots)
library(PerformanceAnalytics) # Used to generate a correlation plot
library(Hmisc) # used to generate histograms for all variables in data frame
library(glmmTMB) # Constructing GLMMs## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.4.1
## Current Matrix version is 1.5.3
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' packagelibrary(MuMIn) # for model selection
library(rphylopic) # add animal silhouettes to graphs
library(broom) # extracting odds ratios in a tidy formatData and summary statistics
Load detection data
Read in saved and cleaned detection data for the 6 LUs from 2021-2022 and 2022-2023 e.g., the 1_ACME_camera_script_9-2-2024.R.
These are two separate files from the two different fiscal years, so they need to be imported and then merged into one data file for plotting. Since both are stored in the data/processed/ folder we can read them both in as a list with purrr.
# detection data
# read in saved and cleaned detection data from the 1_ACME_camera_script_9-2-2024.R
detections <- file.path('data/processed',
c('OSM_ind_det_2021.csv',
'OSM_ind_det_2022.csv')) %>%
map(~.x %>%
read_csv(.) %>%
# ensure the array, site, species, and event_id read in as factors
mutate_if(is.character,
as.factor)) %>%
# give names to each data frame in list
purrr::set_names('dets_2021',
'dets_2022') # R doesn't like when they are just numbers, you can make it work but it's annoying to call the data frame later so I've called them dets_year## Rows: 6696 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): array, site, species, event_id
## dbl (3): month, year, timediff
## dttm (1): datetime
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 14063 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): array, site, species, event_id
## dbl (3): month, year, timediff
## dttm (1): datetime
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# look at data structure
str(detections)## List of 2
## $ dets_2021: tibble [6,696 × 8] (S3: tbl_df/tbl/data.frame)
## ..$ array : Factor w/ 2 levels "LU2","LU3": 1 1 1 1 1 1 1 1 1 1 ...
## ..$ site : Factor w/ 78 levels "LU2_03","LU2_05",..: 1 1 1 1 1 1 1 1 1 1 ...
## ..$ species : Factor w/ 35 levels "Arctic hare",..: 28 28 35 35 35 35 35 35 35 35 ...
## ..$ datetime: POSIXct[1:6696], format: "2021-07-15 11:40:07" "2022-02-06 15:11:06" ...
## ..$ month : num [1:6696] 7 2 7 8 8 8 8 8 8 9 ...
## ..$ year : num [1:6696] 2021 2022 2021 2021 2021 ...
## ..$ timediff: num [1:6696] NA 296839 NA 28519 13230 ...
## ..$ event_id: Factor w/ 6696 levels "E1000000","E1000001",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ dets_2022: tibble [14,063 × 8] (S3: tbl_df/tbl/data.frame)
## ..$ array : Factor w/ 4 levels "LU01","LU13",..: 1 1 1 1 1 1 1 1 1 1 ...
## ..$ site : Factor w/ 155 levels "LU01_06","LU01_10",..: 1 1 1 1 1 1 1 1 1 1 ...
## ..$ species : Factor w/ 39 levels "ATVer","Beaver",..: 31 31 38 38 38 38 38 38 38 38 ...
## ..$ datetime: POSIXct[1:14063], format: "2022-06-17 10:01:52" "2023-09-10 12:51:15" ...
## ..$ month : num [1:14063] 6 9 6 7 7 7 8 8 8 8 ...
## ..$ year : num [1:14063] 2022 2023 2022 2022 2022 ...
## ..$ timediff: num [1:14063] NA 648166 NA 31847 21429 ...
## ..$ event_id: Factor w/ 14063 levels "E0","E1","E10",..: 1 2 5176 6287 7398 8509 9620 10731 11842 12953 ...Merge data
Now we need to merge these two files to make plotting them easier.
They have the same columns so we could just the base R
rbind() function, but in case there are differences in the
columns in the future, let’s use the cleaner
dplyr::bind_rows() function
# Join two years of detection data
detections_merged <- dplyr::bind_rows(detections$dets_2021,
detections$dets_2022)
# check structure of new data
str(detections_merged)## tibble [20,759 × 8] (S3: tbl_df/tbl/data.frame)
## $ array : Factor w/ 6 levels "LU2","LU3","LU01",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ site : Factor w/ 233 levels "LU2_03","LU2_05",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ species : Factor w/ 42 levels "Arctic hare",..: 28 28 35 35 35 35 35 35 35 35 ...
## $ datetime: POSIXct[1:20759], format: "2021-07-15 11:40:07" "2022-02-06 15:11:06" ...
## $ month : num [1:20759] 7 2 7 8 8 8 8 8 8 9 ...
## $ year : num [1:20759] 2021 2022 2021 2021 2021 ...
## $ timediff: num [1:20759] NA 296839 NA 28519 13230 ...
## $ event_id: Factor w/ 20759 levels "E1000000","E1000001",..: 1 2 3 4 5 6 7 8 9 10 ...Data checks
Check to make sure the data looks good after merging before moving on.
Arrays and sites
Let’s check that there are the correct number of levels for arrays and sites, there should be 6 arrays and 233 sites (155 from 2022-2023 and 78 from 2021-2022).
# this will provide a summary of the levels
str(detections_merged)## tibble [20,759 × 8] (S3: tbl_df/tbl/data.frame)
## $ array : Factor w/ 6 levels "LU2","LU3","LU01",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ site : Factor w/ 233 levels "LU2_03","LU2_05",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ species : Factor w/ 42 levels "Arctic hare",..: 28 28 35 35 35 35 35 35 35 35 ...
## $ datetime: POSIXct[1:20759], format: "2021-07-15 11:40:07" "2022-02-06 15:11:06" ...
## $ month : num [1:20759] 7 2 7 8 8 8 8 8 8 9 ...
## $ year : num [1:20759] 2021 2022 2021 2021 2021 ...
## $ timediff: num [1:20759] NA 296839 NA 28519 13230 ...
## $ event_id: Factor w/ 20759 levels "E1000000","E1000001",..: 1 2 3 4 5 6 7 8 9 10 ...# this will list all unique entries (levels) for each variable
levels(detections_merged$array)## [1] "LU2" "LU3" "LU01" "LU13" "LU15" "LU21"levels(detections_merged$site)## [1] "LU2_03" "LU2_05" "LU2_100" "LU2_101" "LU2_106" "LU2_110"
## [7] "LU2_119" "LU2_123" "LU2_13" "LU2_137" "LU2_143" "LU2_15"
## [13] "LU2_153" "LU2_155" "LU2_17" "LU2_18" "LU2_20" "LU2_21"
## [19] "LU2_24" "LU2_25" "LU2_26" "LU2_27" "LU2_29" "LU2_30"
## [25] "LU2_31" "LU2_32" "LU2_33" "LU2_34" "LU2_36" "LU2_38"
## [31] "LU2_39" "LU2_40" "LU2_42" "LU2_44" "LU2_46" "LU2_47"
## [37] "LU2_50" "LU2_51" "LU2_52" "LU2_54" "LU2_56" "LU2_57"
## [43] "LU3_05" "LU3_07" "LU3_10" "LU3_102" "LU3_103" "LU3_111"
## [49] "LU3_126" "LU3_129" "LU3_13" "LU3_131" "LU3_137" "LU3_14"
## [55] "LU3_147" "LU3_15" "LU3_17" "LU3_18" "LU3_19" "LU3_20"
## [61] "LU3_21" "LU3_25" "LU3_27" "LU3_28" "LU3_32" "LU3_36"
## [67] "LU3_38" "LU3_39" "LU3_40" "LU3_41" "LU3_44" "LU3_45"
## [73] "LU3_46" "LU3_48" "LU3_49" "LU3_50" "LU3_51" "LU3_52"
## [79] "LU01_06" "LU01_10" "LU01_11" "LU01_13" "LU01_22" "LU01_25"
## [85] "LU01_27" "LU01_30" "LU01_32" "LU01_36" "LU01_40" "LU01_41"
## [91] "LU01_43" "LU01_44" "LU01_45" "LU01_46" "LU01_47" "LU01_48"
## [97] "LU01_60" "LU01_63" "LU01_64" "LU01_66" "LU01_67" "LU01_70"
## [103] "LU01_71" "LU01_72" "LU01_73" "LU01_74" "LU01_75" "LU01_76"
## [109] "LU01_77" "LU01_78" "LU01_79" "LU01_80" "LU01_82" "LU01_83"
## [115] "LU01_84" "LU01_85" "LU01_86" "LU13_03" "LU13_05" "LU13_06"
## [121] "LU13_08" "LU13_11" "LU13_12" "LU13_128" "LU13_13" "LU13_131"
## [127] "LU13_14" "LU13_15" "LU13_16" "LU13_17" "LU13_18" "LU13_19"
## [133] "LU13_20" "LU13_21" "LU13_22" "LU13_26" "LU13_27" "LU13_30"
## [139] "LU13_32" "LU13_33" "LU13_34" "LU13_35" "LU13_36" "LU13_37"
## [145] "LU13_38" "LU13_41" "LU13_43" "LU13_45" "LU13_47" "LU13_49"
## [151] "LU13_51" "LU13_52" "LU13_53" "LU13_55" "LU13_56" "LU13_57"
## [157] "LU13_59" "LU13_70" "LU15_01" "LU15_02" "LU15_03" "LU15_04"
## [163] "LU15_07" "LU15_08" "LU15_09" "LU15_10" "LU15_11" "LU15_12"
## [169] "LU15_14" "LU15_15" "LU15_16" "LU15_17" "LU15_18" "LU15_19"
## [175] "LU15_20" "LU15_21" "LU15_22" "LU15_23" "LU15_24" "LU15_25"
## [181] "LU15_26" "LU15_27" "LU15_28" "LU15_29" "LU15_30" "LU15_31"
## [187] "LU15_32" "LU15_34" "LU15_36" "LU15_37" "LU15_40" "LU15_41"
## [193] "LU15_43" "LU15_44" "LU15_46" "LU15_58" "LU15_61" "LU21_06"
## [199] "LU21_09" "LU21_10" "LU21_100" "LU21_105" "LU21_106" "LU21_107"
## [205] "LU21_109" "LU21_114" "LU21_116" "LU21_119" "LU21_122" "LU21_126"
## [211] "LU21_14" "LU21_153" "LU21_16" "LU21_164" "LU21_21" "LU21_23"
## [217] "LU21_27" "LU21_32" "LU21_36" "LU21_41" "LU21_52" "LU21_56"
## [223] "LU21_57" "LU21_59" "LU21_63" "LU21_68" "LU21_74" "LU21_78"
## [229] "LU21_82" "LU21_871" "LU21_93" "LU21_97" "LU21_98"Everything looks good.
NAs
Let’s ensure no NAs were introduced where there shouldn’t be during the joining process.
# check for NAs introduced during data merge
summary(detections_merged)## array site species
## LU2 :4711 LU01_47: 405 White-tailed deer:6141
## LU3 :1985 LU01_84: 396 Snowshoe hare :4571
## LU01:6283 LU01_66: 383 Red squirrel :2201
## LU13:1972 LU01_27: 320 Black bear :1997
## LU15:2951 LU2_39 : 310 Coyote :1285
## LU21:2857 LU2_50 : 292 Moose : 693
## (Other):18653 (Other) :3871
## datetime month year
## Min. :2021-07-08 10:21:40.00 Min. : 1.000 Min. :2021
## 1st Qu.:2022-04-25 01:14:32.50 1st Qu.: 5.000 1st Qu.:2022
## Median :2022-11-03 19:10:30.00 Median : 8.000 Median :2022
## Mean :2022-10-09 18:09:20.27 Mean : 7.204 Mean :2022
## 3rd Qu.:2023-05-05 01:40:45.50 3rd Qu.: 9.000 3rd Qu.:2023
## Max. :2023-10-02 11:08:57.00 Max. :12.000 Max. :2023
##
## timediff event_id
## Min. : 30 E1000000: 1
## 1st Qu.: 1427 E1000001: 1
## Median : 5270 E1000002: 1
## Mean : 30292 E1000003: 1
## 3rd Qu.: 18443 E1000004: 1
## Max. :653883 E1000005: 1
## NA's :2446 (Other) :20753THe only NAs are in the timediff column which is what we expect since any of the first observations won’t have a a value for timediff. If you are confused by this re-visit the 1_ACME_camera_script.
Data Formatting
In order to get plots that have the same formatting as last years’ report we have to do a bit of data formatting. First we need to make sure we are including the same relevant species (some were ignored for last years’ report or grouped together).
Last years report had the following species
- white-tailed deer
- snowshoe hare
- black bear
- coyote
- red squirrel
- fisher
- unknown
- moose
- lynx
- spruce grouse
- red fox
- striped skunk
- ruffed grouse
- owl
- grey wolf
- domestic dog
- cougar
- raven
- other
- mule deer
And they grouped all humans except for staff as ‘Humans’. Let’s look at the species we have in the combined years of data and try to format it the same way.
detections_merged %>%
# group by array and species
group_by(species) %>%
summarise(n = n()) %>%
# have R print everything
print(n = nrow(.))## # A tibble: 42 × 2
## species n
## <fct> <int>
## 1 Arctic hare 1
## 2 ATVer 41
## 3 Beaver 4
## 4 Black bear 1997
## 5 Caribou 115
## 6 Cougar 37
## 7 Coyote 1285
## 8 Domestic dog 15
## 9 Elk 1
## 10 Fisher 257
## 11 Grey jay 66
## 12 Grey wolf 224
## 13 Human 16
## 14 Lynx 525
## 15 Marten 220
## 16 Moose 693
## 17 Mule deer 2
## 18 Other 11
## 19 Other birds 219
## 20 Owl 15
## 21 Raven 10
## 22 Red fox 127
## 23 Red squirrel 2201
## 24 Ruffed grouse 90
## 25 Snowmobiler 16
## 26 Snowshoe hare 4571
## 27 Spruce grouse 93
## 28 Staff 461
## 29 Striped skunk 65
## 30 Unknown 663
## 31 Unknown canid 76
## 32 Unknown deer 340
## 33 Unknown mustelid 66
## 34 Unknown ungulate 34
## 35 White-tailed deer 6141
## 36 Canada goose 4
## 37 Hunter 1
## 38 Long-tailed weasel 17
## 39 Otter 7
## 40 Porcupine 5
## 41 Short-tailed weasel 19
## 42 Wolverine 8Hmmm there is one instance of an arctic hare, check that this isn’t meant to be a snowshoe hare and fix later if needed.
Now let’s create a new data frame (tibble) to work with for the OSM figure summaries specifically
I personally would lump all the unknown together and all the birds together but for the sake of consistency with last year’s figures we will remove the same entries and keep the birds separate, let’s create a vector of entries to drop
# create vector of entries to drop for plots
species_drop <- c('Staff',
'Unknown deer',
'Unknown ungulate',
'Unknown canid',
'Unknown mustelid',
'Other birds',
'Arctic hare')
# now we can create the new data frame with some changes consistent w/ choices made for 2021-2022
detections_merged <- detections_merged %>%
# for summarizing, lets lump all the recreational humans into "Humans"
mutate(species = recode_factor(species,
"Snowmobiler" = "Human",
"ATVer" = "Human",
'Hunter' = 'Human')) %>%
# remove species we don't want to plot
filter(!species %in% species_drop)
# look at data
str(detections_merged)## tibble [19,562 × 8] (S3: tbl_df/tbl/data.frame)
## $ array : Factor w/ 6 levels "LU2","LU3","LU01",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ site : Factor w/ 233 levels "LU2_03","LU2_05",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ species : Factor w/ 39 levels "Human","Arctic hare",..: 33 33 33 33 33 33 33 33 33 33 ...
## $ datetime: POSIXct[1:19562], format: "2021-07-18 08:59:30" "2021-08-07 04:21:14" ...
## $ month : num [1:19562] 7 8 8 8 8 8 8 9 9 9 ...
## $ year : num [1:19562] 2021 2021 2021 2021 2021 ...
## $ timediff: num [1:19562] NA 28519 13230 4290 7337 ...
## $ event_id: Factor w/ 20759 levels "E1000000","E1000001",..: 3 4 5 6 7 8 9 10 11 12 ...Subset data
We will also want to subset the data by landscape unit (LU) and generate a new data frame for each LU to use for plotting
I’m not great at writing loops, so let’s see how this shit goes… probably bad but who knows
array_frames <- list()
for (i in unique(detections_merged$array)){
#Subset data based on radius
df <- detections_merged %>%
filter(array == i)
# list of dataframes
array_frames <- c(array_frames, list(df))
}
# inspect one data frame
print(array_frames[[1]]) # this is for LU2## # A tibble: 4,592 × 8
## array site species datetime month year timediff event_id
## <fct> <fct> <fct> <dttm> <dbl> <dbl> <dbl> <fct>
## 1 LU2 LU2_03 White-tailed … 2021-07-18 08:59:30 7 2021 NA E1000002
## 2 LU2 LU2_03 White-tailed … 2021-08-07 04:21:14 8 2021 28519. E1000003
## 3 LU2 LU2_03 White-tailed … 2021-08-16 08:51:21 8 2021 13230. E1000004
## 4 LU2 LU2_03 White-tailed … 2021-08-19 08:21:29 8 2021 4290. E1000005
## 5 LU2 LU2_03 White-tailed … 2021-08-24 10:39:23 8 2021 7337. E1000006
## 6 LU2 LU2_03 White-tailed … 2021-08-26 09:17:02 8 2021 2797. E1000007
## 7 LU2 LU2_03 White-tailed … 2021-08-31 19:13:33 8 2021 7796 E1000008
## 8 LU2 LU2_03 White-tailed … 2021-09-10 05:03:31 9 2021 13550. E1000009
## 9 LU2 LU2_03 White-tailed … 2021-09-16 17:10:41 9 2021 9363. E1000010
## 10 LU2 LU2_03 White-tailed … 2021-09-16 19:19:24 9 2021 127 E1000011
## # ℹ 4,582 more rows… I think this worked
Now let’s change names of list items using purrr, couldn’t figure out how to name them in the loop, you don’t necessarily need to do this because we change the names in the next section, but I like having things named
array_frames <- array_frames %>%
purrr::set_names('LU02',
'LU03',
'LU01',
'LU13',
'LU15',
'LU21')
# inspect each data frame
head(array_frames$LU01)## # A tibble: 6 × 8
## array site species datetime month year timediff event_id
## <fct> <fct> <fct> <dttm> <dbl> <dbl> <dbl> <fct>
## 1 LU01 LU01_06 White-tailed … 2022-06-18 11:09:19 6 2022 NA E2
## 2 LU01 LU01_06 White-tailed … 2022-07-10 13:56:10 7 2022 31847. E3
## 3 LU01 LU01_06 White-tailed … 2022-07-25 11:04:44 7 2022 21429. E4
## 4 LU01 LU01_06 White-tailed … 2022-07-31 06:38:06 7 2022 8356. E5
## 5 LU01 LU01_06 White-tailed … 2022-08-01 09:45:28 8 2022 1627. E6
## 6 LU01 LU01_06 White-tailed … 2022-08-01 15:51:01 8 2022 364. E7Detection plots
Detection data
Now we can apply the same data formatting for each LUs’ data frame using purrr.
We want to count the number of independent detections per species per LU to use in the detection plots
# apply the same formatting to each LU data frame using purrr map
detection_data <- array_frames %>%
purrr::map(
~.x %>%
# group by species
group_by(species) %>%
# calculate a column with unique accounts of each species
mutate(count = n_distinct(event_id)) %>%
# keep just the columns we need
select(species, count) %>%
# keep only unique (distinct) rows so we should be left with one row per species, this helps with plotting later if you don't do it ggplot will try to count and plot each row it's annoying
distinct()) %>%
# set names of list objects
purrr::set_names('Detections LU02',
'Detections LU03',
'Detections LU01',
'Detections LU13',
'Detections LU15',
'Detections LU21')Detection plots
Now to graph independent detections for each LU using purrr, this avoids a TON of code repetition needed to plot each one individually
We use purrr::imap() instead of
purrr::map() because imap maintains the variable names in
our list (e.g. Detections LU01, Detections LU13, etc.) which we can then
use to title each plot.
Within purrr::imap() we just paste the code we would use
for a single ggplot since all the graphical elements (except the title
which we change with the file name [.y]) are the same
# create object detection plots which uses the detection_data list (w/ all 4 LUs)
detection_plots <- detection_data %>%
# use imap instead of map as it allows us to use .y to paste the list element names as the plot titles later
purrr::imap(
~.x %>%
# now just copy and paste the ggplot code for the detection graphs
ggplot(.,
aes(x = reorder(species, count), y = count)) +
# plot as bar graph using geom_col so we don't have to provide a y aesthetic
geom_col() +
# switch the x and y axis
coord_flip() +
# add the number of detections at the end of each bar
geom_text(aes(label = count),
color = "black",
size = 3,
hjust = -0.3,
vjust = 0.2) +
# label x and y axis with informative titles
labs(x = 'Species',
y = 'Number of Independent (30 min) Detections') +
# add title to plot with LU name the .y will take the name of whatever you named each list element in the detection_data list, so make sure this name is what you want on the ggtitle
ggtitle(.y) +
# set the theme
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)))
# view plots, this will print each in it's own window so you have to scroll back in the plot viewer pane to look at each one
detection_plots## $`Detections LU02`
##
## $`Detections LU03`
##
## $`Detections LU01`
##
## $`Detections LU13`
##
## $`Detections LU15`
##
## $`Detections LU21`
Save detection plots
Now we want to save these plots in case we need each individual one (we will combine the detection and naive occ plots into a single figure for each LU later and use those for the OSM report, but we may want these standalone plots later so let’s save them while they are here).
We can save all the plots from the purrr iteration above using
purrr::imap. imap is used instead of map because it allows
us to retain the list object names (plot names) to paste as the file
name with the .y command.
IMPORTANT if you are using this code for a future github repo, DO NOT use .tiff as the file extension. This will cause issues when trying to push any changes to the github repo as the files are too large to meet githubs requirements
# save plots only use if needed
purrr::imap(
detection_plots,
~ggsave(.x,
file = paste0("figures/OSM_",
.y,
'.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
dpi = 600,
width = 11,
height = 9,
units = 'in'))Naive occupancy
Data
We also need to alter the detection data a bit to use for naive occupancy plots.
We will use the individual LU detection data like we did before and
use purrr::map() to apply the dame data formatting to all 4
data frames.
Here we want to calculate the total number of sites in each LU, the number of sites each species was detected at in each LU and then use both those numbers to calculate naive occupancy for each species in each LU
# First we need to alter the data frame a bit for these plots, let's create a data frame for each LU (I couldn't figure out how to do this without assigning individual data frames for each UGH)
# apply the same formatting to each data frame using purrr
occupancy_data <- array_frames %>%
purrr::map(
~.x %>%
# calculate the total number of sites for each LU
mutate(total_sites = n_distinct(site)) %>%
# group by species to calculate the number of sites each spp occurred at
group_by(species) %>%
# add columns to count the number of sites each spp occurred at and then the naive occupancy
reframe(count = n_distinct(site),
naive_occ = count/total_sites,
ind_det = n_distinct(event_id)) %>%
# keep just the columns we need
select(species, naive_occ, ind_det) %>%
# keep only unique (distinct) rows so we should be left with one row per species, this helps with plotting
distinct()) %>%
purrr::set_names('Naive Occupancy LU02',
'Naive Occupancy LU03',
'Naive Occupancy LU01',
'Naive Occupancy LU13',
'Naive Occupancy LU15',
'Naive Occupancy LU21')Occupancy plots
Now we can graph naive occupancy for each LU using purrr, and as with the detection plots this saves a massive amount of coding using purrr to run an iteration on the data files and produce four plots at once instead of copying and pasting code for each individually
# create object occupancy_plots which uses the occupancy_data list (w/ all 4 LUs)
occupancy_plots <- occupancy_data %>%
# use imap instead of map as it allows us to use .y to paste the list element names as the plot titles later
purrr::imap(
~.x %>%
# now just copy and paste the ggplot code for the occupancy graphs
ggplot(.,
aes(x = fct_reorder(species,
ind_det), # this reorders the species so they match the order of the detection plot which makes it better for viewing when the plots are arranged together in 1 figure for each LU
y = naive_occ)) +
# plot as bars using geom_col() which uses stat = 'identity', instead of geom_bar() which will count the rows in each group and plot that instead of naive occ
geom_col() +
# flip x and y axis
coord_flip() +
# add text to end of bars that provides naive occ value
geom_text(aes(label = round(naive_occ, 2)),
size = 3,
hjust = -0.3,
vjust = 0.2) +
# relabel x and y axis and title
labs(x = 'Species',
y = 'Proportion of Sites With At Least One Detection') +
# set plot title using .y (name of list object)
ggtitle(.y) +
# set. theme elements
theme_classic()+
theme(plot.title = element_text(hjust = 0.5)))
# view plots
occupancy_plots## $`Naive Occupancy LU02`
##
## $`Naive Occupancy LU03`
##
## $`Naive Occupancy LU01`
##
## $`Naive Occupancy LU13`
##
## $`Naive Occupancy LU15`
##
## $`Naive Occupancy LU21`
Save occupancy plots
As with the detection plots, we might want these individual plots
later for something so we can use purrr::imap() to save
them to the figures folder
Again avoid using the .tiff extension in github
# save plots
purrr::imap(
occupancy_plots,
~ggsave(.x,
file = paste0("figures/OSM_",
.y,
'.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
dpi = 600,
width = 11,
height = 9,
units = 'in'))Final combined plots for report
The previous year’s report had a figure for each LU with the
detections plot on the top and the occupancy plot on the bottom so we
will recreate these for this year using ggarrange().
Unfortunately I could not figure out how to do this in purrr to reduce coding but luckily it isn’t too much repetition
# not sure I know how to do the following section in purrr just yet, but we've saved a ton of coding so far and it doesn't take much to arrange each of these individually
# LU2
LU02_det_occ_plots <- ggarrange(detection_plots$`Detections LU02`,
occupancy_plots$`Naive Occupancy LU02`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU02_det_occ_plots
# LU3
LU03_det_occ_plots <- ggarrange(detection_plots$`Detections LU03`,
occupancy_plots$`Naive Occupancy LU03`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU03_det_occ_plots
# LU1
# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU01_det_occ_plots <- ggarrange(detection_plots$`Detections LU01`,
occupancy_plots$`Naive Occupancy LU01`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU01_det_occ_plots
# LU13
# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU13_det_occ_plots <- ggarrange(detection_plots$`Detections LU13`,
occupancy_plots$`Naive Occupancy LU13`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU13_det_occ_plots
# LU15
# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU15_det_occ_plots <- ggarrange(detection_plots$`Detections LU15`,
occupancy_plots$`Naive Occupancy LU15`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU15_det_occ_plots
# LU21
# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU21_det_occ_plots <- ggarrange(detection_plots$`Detections LU21`,
occupancy_plots$`Naive Occupancy LU21`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU21_det_occ_plots
We can however, save all the figures again using purrr
# save all figures at once using purrr
final_det_occ_plots <- list(LU02_det_occ_plots,
LU03_det_occ_plots,
LU01_det_occ_plots,
LU13_det_occ_plots,
LU15_det_occ_plots,
LU21_det_occ_plots) %>%
purrr::set_names('LU02_det_occ_plots',
'LU03_det_occ_plots',
'LU01_det_occ_plots',
'LU13_det_occ_plots',
'LU15_det_occ_plots',
'LU21_det_occ_plots') %>%
purrr::imap(
~ggsave(.x,
file = paste0("figures/OSM_",
.y,
'.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
dpi = 600,
width = 12,
height = 15,
units = 'in'))Finish with detection data
Before proceeding let’s clear the objects currently in our environment since we don’t need them for the analysis
rm(list = ls(all.names = TRUE))Global Analysis prep
Now we can start the analysis prep.
First we need to read in the proportional detection (response metrics) and covariate (explanatory metrics) data files for all 6 LUs (fiscal years 2021-2022 and 2022-2023)
Response metrics
Read in data files
We haven’t joined the response metrics from 2021-2022 and 2022-2023 in any previous scripts yet, so let’s read those in using purrr and then join them and take a look at the data to make sure it looks good.
# response metric (proportional detections from the from the ACME_camera_script_9-2-2024.R or .Rmd)
prop_detections <- file.path('data/processed',
c('OSM_proportional_detections_2021.csv',
'OSM_proportional_detections_2022.csv')) %>%
map(~.x %>%
read_csv(.,
# set the column types to read in correctly
col_types = cols(site = col_factor(),
.default = col_number()))) %>%
# give names to each data frame in list
purrr::set_names('prop_dets_2021',
'prop_dets_2022') # R doesn't like when they are just numbers, you can make it work but it's annoying to call the data frame later so I've called them prop_dets_year
# check variable structure
str(prop_detections)## List of 2
## $ prop_dets_2021: spc_tbl_ [78 × 23] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## ..$ site : Factor w/ 78 levels "LU2_03","LU2_05",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..$ black_bear : num [1:78] 2 2 1 3 2 3 4 3 1 2 ...
## ..$ coyote : num [1:78] 1 0 3 5 6 3 8 0 3 4 ...
## ..$ fisher : num [1:78] 2 0 3 1 0 2 1 1 1 0 ...
## ..$ snowshoe_hare : num [1:78] 1 0 2 5 6 9 14 0 7 3 ...
## ..$ white-tailed_deer : num [1:78] 7 3 6 7 6 7 9 6 5 13 ...
## ..$ cougar : num [1:78] 0 1 0 0 0 0 1 0 0 1 ...
## ..$ lynx : num [1:78] 0 0 3 2 0 4 6 0 1 0 ...
## ..$ red_fox : num [1:78] 0 0 3 0 0 0 0 1 0 0 ...
## ..$ moose : num [1:78] 0 0 0 0 1 0 3 3 0 1 ...
## ..$ grey_wolf : num [1:78] 0 0 0 0 0 0 0 0 1 0 ...
## ..$ caribou : num [1:78] 0 0 0 0 0 0 0 0 0 0 ...
## ..$ absent_black_bear : num [1:78] 3 3 4 2 3 8 7 2 4 9 ...
## ..$ absent_coyote : num [1:78] 6 7 4 2 1 11 6 7 4 10 ...
## ..$ absent_fisher : num [1:78] 5 7 4 6 7 12 13 6 6 14 ...
## ..$ absent_snowshoe_hare : num [1:78] 6 7 5 2 1 5 0 7 0 11 ...
## ..$ absent_white-tailed_deer: num [1:78] 0 4 1 0 1 7 5 1 2 1 ...
## ..$ absent_cougar : num [1:78] 7 6 7 7 7 14 13 7 7 13 ...
## ..$ absent_lynx : num [1:78] 7 7 4 5 7 10 8 7 6 14 ...
## ..$ absent_red_fox : num [1:78] 7 7 4 7 7 14 14 6 7 14 ...
## ..$ absent_moose : num [1:78] 7 7 7 7 6 14 11 4 7 13 ...
## ..$ absent_grey_wolf : num [1:78] 7 7 7 7 7 14 14 7 6 14 ...
## ..$ absent_caribou : num [1:78] 7 7 7 7 7 14 14 7 7 14 ...
## ..- attr(*, "spec")=
## .. .. cols(
## .. .. .default = col_number(),
## .. .. site = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
## .. .. black_bear = col_number(),
## .. .. coyote = col_number(),
## .. .. fisher = col_number(),
## .. .. snowshoe_hare = col_number(),
## .. .. `white-tailed_deer` = col_number(),
## .. .. cougar = col_number(),
## .. .. lynx = col_number(),
## .. .. red_fox = col_number(),
## .. .. moose = col_number(),
## .. .. grey_wolf = col_number(),
## .. .. caribou = col_number(),
## .. .. absent_black_bear = col_number(),
## .. .. absent_coyote = col_number(),
## .. .. absent_fisher = col_number(),
## .. .. absent_snowshoe_hare = col_number(),
## .. .. `absent_white-tailed_deer` = col_number(),
## .. .. absent_cougar = col_number(),
## .. .. absent_lynx = col_number(),
## .. .. absent_red_fox = col_number(),
## .. .. absent_moose = col_number(),
## .. .. absent_grey_wolf = col_number(),
## .. .. absent_caribou = col_number()
## .. .. )
## ..- attr(*, "problems")=<externalptr>
## $ prop_dets_2022: spc_tbl_ [154 × 25] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## ..$ site : Factor w/ 154 levels "LU01_06","LU01_10",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..$ black_bear : num [1:154] 7 3 4 7 8 9 4 5 7 7 ...
## ..$ coyote : num [1:154] 4 4 8 10 11 9 11 0 9 4 ...
## ..$ fisher : num [1:154] 4 3 3 3 2 1 1 2 0 3 ...
## ..$ moose : num [1:154] 3 2 5 9 1 0 2 4 1 0 ...
## ..$ snowshoe_hare : num [1:154] 4 1 3 0 8 2 2 0 12 4 ...
## ..$ white-tailed_deer : num [1:154] 12 5 12 12 13 14 15 9 12 10 ...
## ..$ cougar : num [1:154] 0 0 1 0 1 0 0 0 0 0 ...
## ..$ grey_wolf : num [1:154] 0 0 2 0 0 0 1 0 0 0 ...
## ..$ lynx : num [1:154] 0 0 1 0 1 1 0 0 0 2 ...
## ..$ red_fox : num [1:154] 0 0 2 0 0 0 0 0 4 0 ...
## ..$ wolverine : num [1:154] 0 0 0 0 0 0 0 0 0 0 ...
## ..$ caribou : num [1:154] 0 0 0 0 0 0 0 0 0 0 ...
## ..$ absent_black_bear : num [1:154] 5 3 8 5 4 3 8 7 5 5 ...
## ..$ absent_coyote : num [1:154] 10 1 6 5 3 5 4 15 6 11 ...
## ..$ absent_fisher : num [1:154] 10 2 11 12 12 13 14 13 15 12 ...
## ..$ absent_moose : num [1:154] 11 3 9 6 13 14 13 11 14 15 ...
## ..$ absent_snowshoe_hare : num [1:154] 10 4 11 15 6 12 13 15 3 11 ...
## ..$ absent_white-tailed_deer: num [1:154] 2 0 2 3 1 0 0 6 3 5 ...
## ..$ absent_cougar : num [1:154] 14 5 13 15 13 14 15 15 15 15 ...
## ..$ absent_grey_wolf : num [1:154] 14 5 12 15 14 14 14 15 15 15 ...
## ..$ absent_lynx : num [1:154] 14 5 13 15 13 13 15 15 15 13 ...
## ..$ absent_red_fox : num [1:154] 14 5 12 15 14 14 15 15 11 15 ...
## ..$ absent_wolverine : num [1:154] 14 5 14 15 14 14 15 15 15 15 ...
## ..$ absent_caribou : num [1:154] 14 5 14 15 14 14 15 15 15 15 ...
## ..- attr(*, "spec")=
## .. .. cols(
## .. .. .default = col_number(),
## .. .. site = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
## .. .. black_bear = col_number(),
## .. .. coyote = col_number(),
## .. .. fisher = col_number(),
## .. .. moose = col_number(),
## .. .. snowshoe_hare = col_number(),
## .. .. `white-tailed_deer` = col_number(),
## .. .. cougar = col_number(),
## .. .. grey_wolf = col_number(),
## .. .. lynx = col_number(),
## .. .. red_fox = col_number(),
## .. .. wolverine = col_number(),
## .. .. caribou = col_number(),
## .. .. absent_black_bear = col_number(),
## .. .. absent_coyote = col_number(),
## .. .. absent_fisher = col_number(),
## .. .. absent_moose = col_number(),
## .. .. absent_snowshoe_hare = col_number(),
## .. .. `absent_white-tailed_deer` = col_number(),
## .. .. absent_cougar = col_number(),
## .. .. absent_grey_wolf = col_number(),
## .. .. absent_lynx = col_number(),
## .. .. absent_red_fox = col_number(),
## .. .. absent_wolverine = col_number(),
## .. .. absent_caribou = col_number()
## .. .. )
## ..- attr(*, "problems")=<externalptr>Merge data files
Now we need to merge the two data files for analysis. We can do this with dplyr.
# merge the proportional detections files so there are rows for both fiscal years
prop_dets_all <- dplyr::bind_rows(prop_detections$prop_dets_2021,
prop_detections$prop_dets_2022)
print(prop_dets_all)## # A tibble: 232 × 25
## site black_bear coyote fisher snowshoe_hare `white-tailed_deer` cougar lynx
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 LU2_… 2 1 2 1 7 0 0
## 2 LU2_… 2 0 0 0 3 1 0
## 3 LU2_… 1 3 3 2 6 0 3
## 4 LU2_… 3 5 1 5 7 0 2
## 5 LU2_… 2 6 0 6 6 0 0
## 6 LU2_… 3 3 2 9 7 0 4
## 7 LU2_… 4 8 1 14 9 1 6
## 8 LU2_… 3 0 1 0 6 0 0
## 9 LU2_… 1 3 1 7 5 0 1
## 10 LU2_… 2 4 0 3 13 1 0
## # ℹ 222 more rows
## # ℹ 17 more variables: red_fox <dbl>, moose <dbl>, grey_wolf <dbl>,
## # caribou <dbl>, absent_black_bear <dbl>, absent_coyote <dbl>,
## # absent_fisher <dbl>, absent_snowshoe_hare <dbl>,
## # `absent_white-tailed_deer` <dbl>, absent_cougar <dbl>, absent_lynx <dbl>,
## # absent_red_fox <dbl>, absent_moose <dbl>, absent_grey_wolf <dbl>,
## # absent_caribou <dbl>, wolverine <dbl>, absent_wolverine <dbl>Format data
This looks good except since there were no wolverines in the first fiscal year of monitoring (LU02 and LU03) those columns have NAs for both arrays, we want to replace those NAs with Zeros and move the wolverine column to the correct location
Let’s do that now.
prop_dets_all <- prop_dets_all %>%
# replace NAs introduced from joining data to zeros
replace(is.na(.),
0) %>%
# relocate wolverine column
relocate(.,
wolverine,
.after = caribou)
# check data
head(prop_dets_all)## # A tibble: 6 × 25
## site black_bear coyote fisher snowshoe_hare `white-tailed_deer` cougar lynx
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 LU2_03 2 1 2 1 7 0 0
## 2 LU2_05 2 0 0 0 3 1 0
## 3 LU2_1… 1 3 3 2 6 0 3
## 4 LU2_1… 3 5 1 5 7 0 2
## 5 LU2_1… 2 6 0 6 6 0 0
## 6 LU2_1… 3 3 2 9 7 0 4
## # ℹ 17 more variables: red_fox <dbl>, moose <dbl>, grey_wolf <dbl>,
## # caribou <dbl>, wolverine <dbl>, absent_black_bear <dbl>,
## # absent_coyote <dbl>, absent_fisher <dbl>, absent_snowshoe_hare <dbl>,
## # `absent_white-tailed_deer` <dbl>, absent_cougar <dbl>, absent_lynx <dbl>,
## # absent_red_fox <dbl>, absent_moose <dbl>, absent_grey_wolf <dbl>,
## # absent_caribou <dbl>, absent_wolverine <dbl>Looks good!
Save data
Let’s save the merged and formatted detection data from 2021 and 2022 for future use
write_csv(prop_dets_all,
'data/processed/OSM_proportional_detections_merged_2021_2022.csv')Covariates
In the previous script, 2_ACME_landscape_covariate_exploration_script.Rmd we joined the two fiscal years of data and grouped the covariates for analysis and saved this data as a csv file, so we can read in this file now and we shouldn’t have to do any further formatting at the moment
Read data
We will check the data structure after reading in the file just to make sure everything looks good.
covariates_all <- read_csv('data/processed/OSM_covariates_grouped_2021_2022.csv',
# set the column types to read in correctly
col_types = cols(array = col_factor(),
camera = col_factor(),
site = col_factor(),
buff_dist = col_factor(),
.default = col_number()))## Warning: The following named parsers don't match the column names: camerastr(covariates_all)## spc_tbl_ [4,660 × 19] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ array : Factor w/ 6 levels "LU13","LU15",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ site : Factor w/ 233 levels "LU13_18","LU13_15",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ buff_dist : Factor w/ 20 levels "250","500","750",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ harvest : num [1:4660] 0 0 0.687 0.337 0 ...
## $ pipeline : num [1:4660] 0 0.068 0 0 0.0301 ...
## $ roads : num [1:4660] 0 0.0174 0 0 0 ...
## $ seismic_lines : num [1:4660] 0 0.03277 0 0.00889 0.01144 ...
## $ seismic_lines_3D : num [1:4660] 0 0 0 0 0.0523 ...
## $ trails : num [1:4660] 0.00588 0.0028 0 0.01591 0 ...
## $ transmission_lines: num [1:4660] 0.0642 0 0 0 0.091 ...
## $ veg_edges : num [1:4660] 0 0.0858 0 0 0 ...
## $ wells : num [1:4660] 0 0 0 0 0.0322 ...
## $ lc_grassland : num [1:4660] 0.193 0.348 0 0 0.178 ...
## $ lc_coniferous : num [1:4660] 0.456 0.358 0.186 1 0.822 ...
## $ lc_broadleaf : num [1:4660] 0 0 0 0 0 ...
## $ lc_mixed : num [1:4660] 0 0.101 0.255 0 0 ...
## $ lc_developed : num [1:4660] 0 0.0916 0 0 0 ...
## $ lc_shrub : num [1:4660] 0.316 0 0.559 0 0 ...
## $ osm_industrial : num [1:4660] 0.383 0.157 0 0 0 ...
## - attr(*, "spec")=
## .. cols(
## .. .default = col_number(),
## .. array = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
## .. site = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
## .. buff_dist = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
## .. harvest = col_number(),
## .. pipeline = col_number(),
## .. roads = col_number(),
## .. seismic_lines = col_number(),
## .. seismic_lines_3D = col_number(),
## .. trails = col_number(),
## .. transmission_lines = col_number(),
## .. veg_edges = col_number(),
## .. wells = col_number(),
## .. lc_grassland = col_number(),
## .. lc_coniferous = col_number(),
## .. lc_broadleaf = col_number(),
## .. lc_mixed = col_number(),
## .. lc_developed = col_number(),
## .. lc_shrub = col_number(),
## .. osm_industrial = col_number()
## .. )
## - attr(*, "problems")=<externalptr>Everything looks good!
###Subset data by buffer
We do need to subset the data so we have separate data frames for each buffer width to work with in the analysis AND to explore correlations between variables at each buffer width, as these may very with spatial scales
Let’s use a for loop to subset the data, thanks Andrew!
buffer_frames <- list()
for (i in unique(covariates_all$buff_dist)){
print(i)
# Subset data based on radius
df <- covariates_all %>%
filter(buff_dist == i)
# list of dataframes
buffer_frames <-c (buffer_frames, list(df))
}## [1] "250"
## [1] "500"
## [1] "750"
## [1] "1000"
## [1] "1250"
## [1] "1500"
## [1] "1750"
## [1] "2000"
## [1] "2250"
## [1] "2500"
## [1] "2750"
## [1] "3000"
## [1] "3250"
## [1] "3500"
## [1] "3750"
## [1] "4000"
## [1] "4250"
## [1] "4500"
## [1] "4750"
## [1] "5000"# name list objects so we can extract names for plotting
buffer_frames <- buffer_frames %>%
# absurdly long way to do this but for sake of time fuck it
purrr::set_names('250 meter buffer',
'500 meter buffer',
'750 meter buffer',
'1000 meter buffer',
'1250 meter buffer',
'1500 meter buffer',
'1750 meter buffer',
'2000 meter buffer',
'2250 meter buffer',
'2500 meter buffer',
'2750 meter buffer',
'3000 meter buffer',
'3250 meter buffer',
'3500 meter buffer',
'3750 meter buffer',
'4000 meter buffer',
'4250 meter buffer',
'4500 meter buffer',
'4750 meter buffer',
'5000 meter buffer')Now we have a list with data frames for each buffer width which we can work with later.
Add response metric
Now that we have the covariate data formatted we need to add the response metric (monthly proportional presence/absence) to the data frames
osm_final_df_2021_2022 <- buffer_frames %>%
purrr::map(
~.x %>%
left_join(prop_dets_all,
by = 'site'))Finish with data formatting
Let’s remove the objects we no longer need from the environment to keep our work space clean
rm(covariates_all,
prop_detections,
df,
i)Global Analysis
Now we are going to run a global model which includes all HFI and LC variables that at first glance (will do a more thorough check later) seem to have enough data to include as covariates for each buffer width, and then we will compare these models see which buffer width best fit the data for each species. After that we will optimize models so they don’t includes any variables that are highly correlated.
Edit: we added subset analsyis (anthropogenic and landscape) where we subset the data into anthropogenic features (HFI) and landscape (landcover classes) to see how optimum buffer size was influenced based on the type of variables included in the models. For these subset analyses we also looked more closely at which variables were correlated and have adjusted the global models accordingly (i.e. dropped/merged the same variables) to make them more accurate and comparable to the subset models. If there are variables commented out of the models, these were used in the initial run-through of the analysis but were dropped in the subset models so I’ve dropped them accordingly in a re-run of the global models.
This almost means we have to add one quick data formatting step which wasn’t here in the initial run through of the data
osm_final_df_2021_2022 <- osm_final_df_2021_2022 %>%
# use purr so all changes are made to all data frames
purrr::map(
~.x %>%
# mutate data to merge variables that were correlated and similar enough to count as one for the purposes of this analysis (full explanation in subset analysis prep section)
mutate(pipeline_transmission_lines = pipeline + transmission_lines,
lc_forest = lc_coniferous + lc_broadleaf + lc_mixed))
# view structure of one data frame
str(osm_final_df_2021_2022$`250 meter buffer`)## tibble [233 × 45] (S3: tbl_df/tbl/data.frame)
## $ array : Factor w/ 6 levels "LU13","LU15",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ site : Factor w/ 233 levels "LU13_18","LU13_15",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ buff_dist : Factor w/ 20 levels "250","500","750",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ harvest : num [1:233] 0 0 0.687 0.337 0 ...
## $ pipeline : num [1:233] 0 0.068 0 0 0.0301 ...
## $ roads : num [1:233] 0 0.0174 0 0 0 ...
## $ seismic_lines : num [1:233] 0 0.03277 0 0.00889 0.01144 ...
## $ seismic_lines_3D : num [1:233] 0 0 0 0 0.0523 ...
## $ trails : num [1:233] 0.00588 0.0028 0 0.01591 0 ...
## $ transmission_lines : num [1:233] 0.0642 0 0 0 0.091 ...
## $ veg_edges : num [1:233] 0 0.0858 0 0 0 ...
## $ wells : num [1:233] 0 0 0 0 0.0322 ...
## $ lc_grassland : num [1:233] 0.193 0.348 0 0 0.178 ...
## $ lc_coniferous : num [1:233] 0.456 0.358 0.186 1 0.822 ...
## $ lc_broadleaf : num [1:233] 0 0 0 0 0 ...
## $ lc_mixed : num [1:233] 0 0.101 0.255 0 0 ...
## $ lc_developed : num [1:233] 0 0.0916 0 0 0 ...
## $ lc_shrub : num [1:233] 0.316 0 0.559 0 0 ...
## $ osm_industrial : num [1:233] 0.383 0.157 0 0 0 ...
## $ black_bear : num [1:233] 1 4 6 2 0 3 1 2 1 1 ...
## $ coyote : num [1:233] 9 8 5 3 2 7 9 0 0 1 ...
## $ fisher : num [1:233] 0 0 0 1 0 0 0 0 0 0 ...
## $ snowshoe_hare : num [1:233] 11 9 0 0 2 3 8 3 0 3 ...
## $ white-tailed_deer : num [1:233] 4 5 3 0 0 3 4 0 0 0 ...
## $ cougar : num [1:233] 0 0 0 0 0 0 0 0 0 0 ...
## $ lynx : num [1:233] 4 1 0 0 1 4 1 1 0 2 ...
## $ red_fox : num [1:233] 0 0 0 0 0 0 0 0 0 0 ...
## $ moose : num [1:233] 0 0 4 0 1 1 0 1 0 3 ...
## $ grey_wolf : num [1:233] 0 3 0 3 0 0 0 0 0 0 ...
## $ caribou : num [1:233] 0 0 0 0 0 0 0 0 0 0 ...
## $ wolverine : num [1:233] 0 0 0 0 0 0 0 0 0 0 ...
## $ absent_black_bear : num [1:233] 8 5 3 7 9 6 8 7 8 8 ...
## $ absent_coyote : num [1:233] 3 4 7 9 10 5 3 12 12 11 ...
## $ absent_fisher : num [1:233] 12 12 12 11 12 12 12 12 12 12 ...
## $ absent_snowshoe_hare : num [1:233] 1 3 12 12 10 9 4 9 12 9 ...
## $ absent_white-tailed_deer : num [1:233] 8 7 9 12 12 9 8 12 12 12 ...
## $ absent_cougar : num [1:233] 12 12 12 12 12 12 12 12 12 12 ...
## $ absent_lynx : num [1:233] 8 11 12 12 11 8 11 11 12 10 ...
## $ absent_red_fox : num [1:233] 12 12 12 12 12 12 12 12 12 12 ...
## $ absent_moose : num [1:233] 12 12 8 12 11 11 12 11 12 9 ...
## $ absent_grey_wolf : num [1:233] 12 9 12 9 12 12 12 12 12 12 ...
## $ absent_caribou : num [1:233] 12 12 12 12 12 12 12 12 12 12 ...
## $ absent_wolverine : num [1:233] 12 12 12 12 12 12 12 12 12 12 ...
## $ pipeline_transmission_lines: num [1:233] 0.0642 0.068 0 0 0.1211 ...
## $ lc_forest : num [1:233] 0.456 0.459 0.441 1 0.822 ...We don’t need to do ALL the species since many don’t have enough data.
Refer to the 1_ACME_camera_script_9-2-2024.html or .Rmd the plot for proportional monthly detections should provide info on which species we have enough data for, can be found under Response metrics/3.Proportion monthly detections
A brief look at this fig indicates that we have enough for all the mammals in the prop_detections data frame except
- cougar
- wolverine
- caribou??? (may have enough, may not)
Black bear
there is probably a way to shorten the following code to select particular species, I saw Andrew’s for loop in the draft script he wrote but couldn’t quite figure out how to adapt it to my purposes with the data formatted the way I have it, so I did this instead, maybe we can merge approaches later to clean this up if deemed necessary? But it certainly functions for now and is understandable… I think.
Global models
Let’s start with bears and use purrr to create a global model for every buffer distance
Recall purrr::map() is magical for iterations and will
apply all the functions within the map() function to each
item of the list supplied before the the map()
function.
# create models for black bears at each buffer size
black_bear_mods <- osm_final_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
# HFI
scale(harvest) +
# scale(pipeline) +
# scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
# scale(transmission_lines) +
# scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# VEG covariates in numerical order
scale(lc_grassland) +
# scale(lc_coniferous) +
# scale(lc_broadleaf) +
# scale(lc_mixed) +
scale(lc_developed) +
# scale(lc_shrub) +
scale(lc_forest) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))Model selection
We will use the model.sel() function from the
MuMIn package to compare the global models for each buffer
width and see which buffer fits the black bear data best
model.sel(black_bear_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_dvl))
## 250 meter buffer -0.5980 + 0.007994 -0.18820
## 3750 meter buffer -0.5985 + 0.063740 -0.12240
## 4000 meter buffer -0.5989 + 0.066910 -0.11740
## 3500 meter buffer -0.5983 + 0.055590 -0.11950
## 3250 meter buffer -0.5980 + 0.057680 -0.11760
## 4250 meter buffer -0.5986 + 0.074160 -0.11000
## 2750 meter buffer -0.5965 + 0.051590 -0.12580
## 3000 meter buffer -0.5970 + 0.054350 -0.11950
## 2500 meter buffer -0.5953 + 0.043640 -0.12320
## 4500 meter buffer -0.5981 + 0.084460 -0.09583
## 2250 meter buffer -0.5949 + 0.034920 -0.12630
## 2000 meter buffer -0.5958 + 0.035770 -0.12270
## 1750 meter buffer -0.5966 + 0.037730 -0.12030
## 1500 meter buffer -0.5968 + 0.038050 -0.13840
## 4750 meter buffer -0.5983 + 0.084940 -0.08955
## 5000 meter buffer -0.5978 + 0.089220 -0.09663
## 1250 meter buffer -0.5962 + 0.028780 -0.13590
## 1000 meter buffer -0.5953 + 0.019340 -0.10210
## 500 meter buffer -0.5968 + 0.032870 -0.11200
## 750 meter buffer -0.5945 + 0.031270 -0.10860
## cnd(scl(lc_frs)) cnd(scl(lc_grs)) cnd(scl(osm_ind))
## 250 meter buffer -0.117300 -0.001634 0.033050
## 3750 meter buffer -0.018580 -0.033640 -0.110600
## 4000 meter buffer -0.019160 -0.021770 -0.110900
## 3500 meter buffer -0.004563 -0.027070 -0.093470
## 3250 meter buffer 0.005321 -0.027890 -0.080990
## 4250 meter buffer -0.030800 -0.017550 -0.110600
## 2750 meter buffer 0.024430 -0.039480 -0.049160
## 3000 meter buffer 0.016130 -0.035470 -0.064570
## 2500 meter buffer 0.028730 -0.060960 -0.038810
## 4500 meter buffer -0.041340 -0.022060 -0.117700
## 2250 meter buffer 0.039870 -0.065080 -0.024960
## 2000 meter buffer 0.041850 -0.051020 -0.009152
## 1750 meter buffer 0.031950 -0.038190 0.002659
## 1500 meter buffer 0.016520 -0.033480 0.017320
## 4750 meter buffer -0.047390 -0.027060 -0.117000
## 5000 meter buffer -0.058540 -0.030800 -0.116400
## 1250 meter buffer 0.015880 -0.027750 0.034920
## 1000 meter buffer 0.036220 -0.016250 0.013740
## 500 meter buffer -0.022960 0.051660 0.018230
## 750 meter buffer 0.028620 0.007247 0.030200
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 250 meter buffer -0.07255 -0.1169 -0.01783
## 3750 meter buffer -0.11560 -0.1217 -0.12120
## 4000 meter buffer -0.12730 -0.1160 -0.12740
## 3500 meter buffer -0.12510 -0.1137 -0.10580
## 3250 meter buffer -0.13400 -0.1098 -0.09997
## 4250 meter buffer -0.11110 -0.1191 -0.14160
## 2750 meter buffer -0.13670 -0.1125 -0.09440
## 3000 meter buffer -0.13650 -0.1104 -0.09114
## 2500 meter buffer -0.11760 -0.1118 -0.09321
## 4500 meter buffer -0.10940 -0.1254 -0.15660
## 2250 meter buffer -0.10900 -0.1123 -0.09143
## 2000 meter buffer -0.10620 -0.1147 -0.10760
## 1750 meter buffer -0.10940 -0.1199 -0.12240
## 1500 meter buffer -0.09957 -0.1287 -0.10070
## 4750 meter buffer -0.10330 -0.1283 -0.15910
## 5000 meter buffer -0.09421 -0.1359 -0.16130
## 1250 meter buffer -0.08952 -0.1341 -0.06911
## 1000 meter buffer -0.09475 -0.1301 -0.05778
## 500 meter buffer -0.12200 -0.1184 -0.02209
## 750 meter buffer -0.08703 -0.1315 -0.03521
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 250 meter buffer 0.12650 -0.01118 12 -448.961 923.3 0.00 0.678
## 3750 meter buffer 0.10380 0.28510 12 -451.536 928.5 5.15 0.052
## 4000 meter buffer 0.10410 0.28820 12 -451.614 928.7 5.31 0.048
## 3500 meter buffer 0.09673 0.26180 12 -452.141 929.7 6.36 0.028
## 3250 meter buffer 0.09012 0.26160 12 -452.235 929.9 6.55 0.026
## 4250 meter buffer 0.10150 0.27290 12 -452.471 930.4 7.02 0.020
## 2750 meter buffer 0.07929 0.25930 12 -452.474 930.4 7.03 0.020
## 3000 meter buffer 0.07835 0.25970 12 -452.688 930.8 7.45 0.016
## 2500 meter buffer 0.08615 0.24080 12 -452.810 931.0 7.70 0.014
## 4500 meter buffer 0.10040 0.28170 12 -452.817 931.1 7.71 0.014
## 2250 meter buffer 0.09190 0.23010 12 -452.861 931.1 7.80 0.014
## 2000 meter buffer 0.08802 0.19780 12 -452.943 931.3 7.96 0.013
## 1750 meter buffer 0.07880 0.17750 12 -453.054 931.5 8.19 0.011
## 1500 meter buffer 0.07856 0.16480 12 -453.190 931.8 8.46 0.010
## 4750 meter buffer 0.09694 0.27290 12 -453.241 931.9 8.56 0.009
## 5000 meter buffer 0.09090 0.27590 12 -453.402 932.2 8.88 0.008
## 1250 meter buffer 0.08315 0.13020 12 -453.454 932.3 8.99 0.008
## 1000 meter buffer 0.08348 0.11150 12 -453.876 933.2 9.83 0.005
## 500 meter buffer 0.07082 0.08991 12 -454.054 933.5 10.19 0.004
## 750 meter buffer 0.04348 0.12260 12 -455.171 935.8 12.42 0.001
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Looks like the smallest buffer (250) fits the data best for black bears.
Note about re-run global models: when we removed the correlated variables the top model satyed the same (250m) but the weight increased (from 0.159 to 0.678)
Let’s look at this model closer
Model summary
Since we aren’t interpreting the magnitude and direction of effect for individual variables in the model as you normally would, the model summary serves primarily to see that there are no issues with the top model (convergence issues, large SE, etc.) that could cause us to quetion the results of the model.sel() function.
summary(black_bear_mods$`250 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(black_bear, absent_black_bear) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## scale(lc_grassland) + scale(lc_developed) + scale(lc_forest) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 921.9 963.3 -449.0 897.9 220
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.08116 0.2849
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.598017 0.125793 -4.754 1.99e-06 ***
## scale(harvest) 0.007994 0.049036 0.163 0.87050
## scale(seismic_lines) -0.017830 0.052626 -0.339 0.73476
## scale(seismic_lines_3D) -0.116896 0.058268 -2.006 0.04484 *
## scale(trails) 0.126499 0.046332 2.730 0.00633 **
## scale(wells) -0.011184 0.053722 -0.208 0.83509
## scale(osm_industrial) 0.033050 0.051185 0.646 0.51848
## scale(pipeline_transmission_lines) -0.072550 0.063235 -1.147 0.25125
## scale(lc_grassland) -0.001634 0.059234 -0.028 0.97799
## scale(lc_developed) -0.188178 0.065107 -2.890 0.00385 **
## scale(lc_forest) -0.117330 0.070534 -1.663 0.09622 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nothing looks fishy in the model summary for now, we will look at this more closely once we have a true top model.
Code to remove a model
At one point the 250 meter buffer was giving us issues so we had to remove it. After re-extracting the data and re-doing some data formatting this is no longer an issue but I’ve saved the code here in case it’s needed in the future
# create models for black bears at each buffer size
black_bear_mods_no250 <- osm_final_df_2021_2022 %>%
# remove 250 meter buffer width
purrr::discard_at('250 meter buffer') %>%
# use purrr map to fun the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
# HFI
scale(harvest) +
scale(pipeline) +
scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(transmission_lines) +
scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
# VEG covariates in numerical order
scale(lc_grassland) +
scale(lc_coniferous) +
scale(lc_broadleaf) +
scale(lc_mixed) +
scale(lc_developed) +
scale(lc_shrub) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))Subset models
NOTE this code was not used for publication but was an initial idea for how we could further use the data. We only explored subset models for Bears thus far, but have left the headings for the other species in case we want to add in these subset models later
Autocorrelation 250m
Before we can develop model subsets we need to see what variables can be included in the same model at this buffer width.
Let’s use the chart.Correlation() function in the
Performance Analytics package to look at this.
buffer_frames$`250 meter buffer` %>%
select_if(is.numeric) %>%
# use chart.correlation
chart.Correlation(.,
histogram = TRUE,
method = "pearson")
mtext('250 meter buffer', side = 3, line = 3)
> You can click on this fig to zoom in!
List of correlated variables:
- pipeline & transmission_lines 0.53
- roads & veg_edges 0.71
- roads & lc_developed 0.57
Model list
Let’s create another global model without these correlated variables. I’m going to select transmission_lines over pipeline because the summary from earlier showed transmission lines had larger effect on black bear presence, and I’m going to choose to keep roads instead of veg edges and the developed landcover class because we are interested in the effect of roads more than these other two variables.
# global model w/ non-correlated variables
bear_global_250 <- glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
# HFI
scale(harvest) +
scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(transmission_lines) +
scale(wells) +
scale(osm_industrial) +
# VEG covariates in numerical order
scale(lc_grassland) +
scale(lc_coniferous) +
scale(lc_broadleaf) +
scale(lc_mixed) +
scale(lc_shrub) +
# Random effect of array
(1|array),
data = osm_final_df_2021_2022$`250 meter buffer`,
family = 'binomial')
# null model to compare
bear_null_250 <- glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~ 1 +
# Random effect of array
(1|array),
data = osm_final_df_2021_2022$`250 meter buffer`,
family = 'binomial')
# second model is based on linear features providing easier movement through boreal forest
bear_linear_250 <- glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
# HFI
scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(transmission_lines) +
# Random effect of array
(1|array),
data = osm_final_df_2021_2022$`250 meter buffer`,
family = 'binomial')Model selection
Now let’s look at a model selection table with our subset models.
model.sel(bear_global_250,
bear_null_250,
bear_linear_250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## bear_linear_250 -0.5986 +
## bear_global_250 -0.5995 + 0.01223 0.2326
## bear_null_250 -0.5884 +
## cnd(scl(lc_cnf)) cnd(scl(lc_grs)) cnd(scl(lc_mxd))
## bear_linear_250
## bear_global_250 0.1823 0.1248 0.1455
## bear_null_250
## cnd(scl(lc_shr)) cnd(scl(osm_ind)) cnd(scl(rds))
## bear_linear_250 -0.1586
## bear_global_250 0.1891 0.04922 -0.1215
## bear_null_250
## cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns)) cnd(scl(trl))
## bear_linear_250 -0.1311 -0.04436 0.10250
## bear_global_250 -0.1135 -0.04271 0.09811
## bear_null_250
## cnd(scl(trn_lns)) cnd(scl(wll)) df logLik AICc delta weight
## bear_linear_250 -0.09901 7 -449.674 913.8 0.00 0.972
## bear_global_250 -0.10520 -0.01348 15 -444.383 921.0 7.14 0.027
## bear_null_250 2 -463.090 930.2 16.38 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Looks like the linear feature model is best, SO FAR
Top model
250 meter buffer - linear features
Summary
summary(bear_linear_250)## Family: binomial ( logit )
## Formula:
## cbind(black_bear, absent_black_bear) ~ scale(roads) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(transmission_lines) +
## (1 | array)
## Data: osm_final_df_2021_2022$`250 meter buffer`
##
## AIC BIC logLik deviance df.resid
## 913.3 937.5 -449.7 899.3 225
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.05426 0.2329
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.59856 0.10644 -5.624 1.87e-08 ***
## scale(roads) -0.15864 0.05694 -2.786 0.00534 **
## scale(seismic_lines) -0.04436 0.05006 -0.886 0.37556
## scale(seismic_lines_3D) -0.13114 0.05635 -2.327 0.01995 *
## scale(trails) 0.10247 0.04641 2.208 0.02724 *
## scale(transmission_lines) -0.09901 0.05653 -1.751 0.07987 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Odds ratios
Let’s extract the odds ratios for the top model so we can plot them for data vis later.
bear_model_odds <- bear_linear_250 %>%
# extract the coefficients and upper and lower CI
confint() %>%
# format resulting object as a tibble data frame
as_tibble() %>%
# add a column where we can put the feature names
rowid_to_column() %>%
# rename the columns for plotting
rename('lower' = `2.5 %`,
'upper' = `97.5 %`,
'estimate' = Estimate,
'feature' = rowid) %>%
# rename the entries to features, need to look at the order the features are in from the model summary and ensure it matches
mutate(feature = as.factor(feature),
feature = recode(feature,
'1' = 'intercept',
'2' = 'roads',
'3' = 'seismic_lines',
'4' = 'seismic_lines_3D',
'5' = 'trails',
'6' = 'transmission_lines',
'7' = 'intercept_array')) %>%
# remove intercepts
filter(!grepl('intercept',
feature))Plot odds ratios
First let’s get a silhouette for this graphy from phylopic
black_bear_img <- get_phylopic(get_uuid(name = 'Ursus americanus'))Now let’s use ggplot to plot the odds ratios for each feature in the top model
# provide data and mapping aesthetics
ggplot(bear_model_odds, aes(x = feature,
y = estimate)) +
geom_errorbar(aes(ymin = lower,
ymax = upper),
width = 0.4,
linewidth = 0.5,
position = position_dodge(width = 0.9)) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(y = 'Odds ratio') +
scale_x_discrete(labels = c('Roads',
'Seismic Lines',
'3D Seismic Lines',
'Trails',
'Transmission Lines')) +
add_phylopic(black_bear_img,
x = 5.3,
y = 0.2,
ysize = 0.05) +
theme_classic() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90,
hjust = 1),
axis.title = element_text(size = 16),
axis.text = element_text(size = 14))
Let’s repeat this process for each species that we have enough data for.
Caribou
We may or may not have enough data for caribou but let’s try it at least for this preliminary report
We can use the same code from black bears (above) to run global models for each buffer width
And in the same chunk to save time let’s also run the
model.sel() function
Global models
caribou_mods <- osm_final_df_2021_2022 %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(caribou, absent_caribou) ~
# HFI
scale(harvest) +
# scale(pipeline) +
# scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
# scale(transmission_lines) +
# scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# VEG covariates in numerical order
scale(lc_grassland) +
# scale(lc_coniferous) +
# scale(lc_broadleaf) +
# scale(lc_mixed) +
scale(lc_developed) +
# scale(lc_shrub) +
scale(lc_forest) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')# model selection
model.sel(caribou_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_dvl))
## 1000 meter buffer -5.542 + -0.65730 0.06848
## 1250 meter buffer -5.389 + -0.36120 -0.09234
## 750 meter buffer -5.766 + -1.56400 -0.04326
## 1500 meter buffer -5.303 + -0.19770 -0.30870
## 2500 meter buffer -5.604 + -0.16270 0.32760
## 2000 meter buffer -5.379 + -0.21430 0.10330
## 1750 meter buffer -5.283 + -0.21810 -0.31570
## 2750 meter buffer -5.573 + -0.13430 0.30490
## 2250 meter buffer -5.475 + -0.19150 0.24040
## 500 meter buffer -5.750 + -2.30500 -0.05809
## 3000 meter buffer -5.603 + -0.04354 0.44240
## 250 meter buffer -204.200 + -599.00000 -0.30360
## 3250 meter buffer -5.702 + 0.02384 0.61270
## 3500 meter buffer -5.674 + -0.01945 0.62170
## 3750 meter buffer -5.640 + -0.08866 0.62130
## 4750 meter buffer -5.193 + -1.17600 -0.59520
## 5000 meter buffer -5.222 + -1.38100 -0.65940
## 4000 meter buffer -5.528 + -0.19870 0.51390
## 4500 meter buffer -5.133 + -0.89510 -0.25990
## 4250 meter buffer -5.196 + -0.54080 0.02969
## cnd(scl(lc_frs)) cnd(scl(lc_grs)) cnd(scl(osm_ind))
## 1000 meter buffer 0.205300 -0.06587 -1.29300
## 1250 meter buffer 0.190900 -0.15000 -1.14200
## 750 meter buffer 0.253700 -0.06970 -1.02300
## 1500 meter buffer 0.207700 -0.24010 -0.91560
## 2500 meter buffer 0.176600 -0.34140 -1.48800
## 2000 meter buffer 0.192200 -0.27850 -1.50300
## 1750 meter buffer 0.180900 -0.26740 -1.09100
## 2750 meter buffer 0.140400 -0.31460 -1.10700
## 2250 meter buffer 0.192000 -0.34030 -1.62400
## 500 meter buffer 0.159400 -0.13690 -0.31640
## 3000 meter buffer 0.098770 -0.25210 -0.97030
## 250 meter buffer 0.177200 -0.14610 -0.19130
## 3250 meter buffer 0.062340 -0.24330 -0.92670
## 3500 meter buffer 0.046120 -0.26300 -0.81730
## 3750 meter buffer 0.033180 -0.28880 -0.75800
## 4750 meter buffer -0.072110 -0.38770 0.07747
## 5000 meter buffer -0.094550 -0.38140 0.09164
## 4000 meter buffer 0.006902 -0.31770 -0.62650
## 4500 meter buffer -0.032880 -0.36770 -0.07555
## 4250 meter buffer -0.015900 -0.33420 -0.28680
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 1000 meter buffer -0.3002 -0.1386 0.14380
## 1250 meter buffer -0.2207 -0.1474 0.17630
## 750 meter buffer -0.2925 -0.1413 0.21950
## 1500 meter buffer -0.2618 -0.2002 0.08685
## 2500 meter buffer -0.7149 -0.1131 0.16400
## 2000 meter buffer -0.4076 -0.1097 0.10420
## 1750 meter buffer -0.2710 -0.1938 0.02205
## 2750 meter buffer -0.8976 -0.1278 0.11550
## 2250 meter buffer -0.4376 -0.1163 0.17620
## 500 meter buffer -0.1657 -0.1700 0.20930
## 3000 meter buffer -0.9864 -0.1277 0.00994
## 250 meter buffer -0.2073 -0.2941 0.22930
## 3250 meter buffer -0.9571 -0.1512 -0.07157
## 3500 meter buffer -0.9677 -0.1274 -0.08802
## 3750 meter buffer -0.9201 -0.1117 -0.11680
## 4750 meter buffer -0.6927 -0.1673 0.05126
## 5000 meter buffer -0.6612 -0.1574 0.07862
## 4000 meter buffer -0.8016 -0.1044 -0.12110
## 4500 meter buffer -0.7023 -0.1233 -0.01153
## 4250 meter buffer -0.7415 -0.1173 -0.06533
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 1000 meter buffer 0.04559 0.7343 12 -131.820 289.1 0.00 0.796
## 1250 meter buffer 0.05589 0.7962 12 -134.133 293.7 4.63 0.079
## 750 meter buffer 0.03605 0.6896 12 -134.229 293.9 4.82 0.072
## 1500 meter buffer -0.04957 0.8723 12 -135.523 296.5 7.41 0.020
## 2500 meter buffer -0.19600 1.5630 12 -136.308 298.0 8.98 0.009
## 2000 meter buffer -0.06962 1.1630 12 -136.622 298.7 9.60 0.007
## 1750 meter buffer -0.04380 0.9260 12 -136.953 299.3 10.27 0.005
## 2750 meter buffer -0.26150 1.5510 12 -136.958 299.3 10.28 0.005
## 2250 meter buffer -0.09751 1.3500 12 -136.987 299.4 10.33 0.005
## 500 meter buffer -0.04656 0.4975 12 -137.794 301.0 11.95 0.002
## 3000 meter buffer -0.37210 1.4900 12 -138.536 302.5 13.43 0.001
## 250 meter buffer -0.29450 0.3696 12 -138.888 303.2 14.14 0.001
## 3250 meter buffer -0.46890 1.4380 12 -140.040 305.5 16.44 0.000
## 3500 meter buffer -0.50310 1.4150 12 -140.554 306.5 17.47 0.000
## 3750 meter buffer -0.48830 1.3860 12 -141.121 307.7 18.60 0.000
## 4750 meter buffer -0.50800 1.1280 12 -141.749 308.9 19.86 0.000
## 5000 meter buffer -0.52940 1.1220 12 -141.769 309.0 19.90 0.000
## 4000 meter buffer -0.42300 1.3090 12 -142.301 310.0 20.96 0.000
## 4500 meter buffer -0.43140 1.1030 12 -142.883 311.2 22.13 0.000
## 4250 meter buffer -0.41310 1.1760 12 -142.993 311.4 22.35 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)We get a warning that there are some model convergence problems, I expect this is because we don’t have enough data for caribou, so we may choose to remove the caribou results from the manuscript or provide that caveat.
For now let’s examine the top buffer width and see if the model seems reasonable.
Note about re-run global models: when we re-ran the global models with dropped correlated variables our top model for caribou does change (3000m to 1000m) which now matches the top model for anthropogenic features, this could be because there are more anthro features in the global models and/or because of lack of data/model convergence issues. If we see this same trend with other global models it will be worth mentioning as a caveat in the discussion.
Model summary 1000m
Let’s take a closer look at the top model summary
summary(caribou_mods$`1000 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(caribou, absent_caribou) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## scale(lc_grassland) + scale(lc_developed) + scale(lc_forest) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 287.6 329.0 -131.8 263.6 220
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 4.767 2.183
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.54169 1.07006 -5.179 2.23e-07 ***
## scale(harvest) -0.65725 0.64637 -1.017 0.30924
## scale(seismic_lines) 0.14381 0.16616 0.865 0.38678
## scale(seismic_lines_3D) -0.13865 0.21047 -0.659 0.51006
## scale(trails) 0.04559 0.20832 0.219 0.82676
## scale(wells) 0.73427 0.17137 4.285 1.83e-05 ***
## scale(osm_industrial) -1.29295 0.47653 -2.713 0.00666 **
## scale(pipeline_transmission_lines) -0.30020 0.30963 -0.970 0.33227
## scale(lc_grassland) -0.06587 0.16616 -0.396 0.69178
## scale(lc_developed) 0.06848 0.40830 0.168 0.86681
## scale(lc_forest) 0.20534 0.14382 1.428 0.15336
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There’s nothing that catches my eye immediately as being sus about this particular model so it may not have been one with convergence issues. We will keep it in report for now
Subset Models
Add later if deemed necessary.
Coyote
Global models
coyote_mods <- osm_final_df_2021_2022 %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(coyote, absent_coyote) ~
# HFI
scale(harvest) +
# scale(pipeline) +
# scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
# scale(transmission_lines) +
# scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# VEG covariates in numerical order
scale(lc_grassland) +
# scale(lc_coniferous) +
# scale(lc_broadleaf) +
# scale(lc_mixed) +
scale(lc_developed) +
# scale(lc_shrub) +
scale(lc_forest) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(coyote_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_dvl))
## 5000 meter buffer -1.387 + -0.084680 0.3724
## 4750 meter buffer -1.386 + -0.080090 0.3626
## 4500 meter buffer -1.385 + -0.070060 0.3539
## 4250 meter buffer -1.383 + -0.061960 0.3338
## 1250 meter buffer -1.388 + 0.041420 0.3791
## 1500 meter buffer -1.383 + 0.040670 0.3472
## 2250 meter buffer -1.378 + -0.011870 0.3432
## 1750 meter buffer -1.379 + 0.026290 0.3454
## 2000 meter buffer -1.377 + 0.003104 0.3411
## 3000 meter buffer -1.380 + -0.024970 0.3105
## 2500 meter buffer -1.379 + -0.016180 0.3310
## 2750 meter buffer -1.379 + -0.021550 0.3161
## 4000 meter buffer -1.382 + -0.059530 0.3264
## 1000 meter buffer -1.391 + 0.042800 0.3792
## 3250 meter buffer -1.380 + -0.031490 0.3142
## 3500 meter buffer -1.379 + -0.043320 0.3315
## 3750 meter buffer -1.379 + -0.049730 0.3305
## 750 meter buffer -1.393 + 0.073420 0.3437
## 500 meter buffer -1.390 + 0.038750 0.2173
## 250 meter buffer -1.386 + -0.056700 0.1532
## cnd(scl(lc_frs)) cnd(scl(lc_grs)) cnd(scl(osm_ind))
## 5000 meter buffer -0.1600 -0.201600 0.3537
## 4750 meter buffer -0.1688 -0.197700 0.3539
## 4500 meter buffer -0.1750 -0.192300 0.3471
## 4250 meter buffer -0.1808 -0.171200 0.3403
## 1250 meter buffer -0.1454 -0.022650 0.2821
## 1500 meter buffer -0.1459 -0.033380 0.3036
## 2250 meter buffer -0.1491 -0.092860 0.3043
## 1750 meter buffer -0.1416 -0.053430 0.3033
## 2000 meter buffer -0.1376 -0.073700 0.3109
## 3000 meter buffer -0.1720 -0.060970 0.3341
## 2500 meter buffer -0.1560 -0.077810 0.3135
## 2750 meter buffer -0.1623 -0.060870 0.3361
## 4000 meter buffer -0.1865 -0.140900 0.3178
## 1000 meter buffer -0.1430 -0.011130 0.2735
## 3250 meter buffer -0.1786 -0.076440 0.3277
## 3500 meter buffer -0.1844 -0.107400 0.3027
## 3750 meter buffer -0.1920 -0.127300 0.2961
## 750 meter buffer -0.1185 0.007329 0.2878
## 500 meter buffer -0.1127 0.020850 0.2932
## 250 meter buffer -0.1637 -0.022140 0.2168
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 5000 meter buffer 0.1137000 -0.07121 0.2669
## 4750 meter buffer 0.1076000 -0.07461 0.2591
## 4500 meter buffer 0.0934800 -0.07683 0.2565
## 4250 meter buffer 0.0789600 -0.07192 0.2365
## 1250 meter buffer 0.0194300 -0.14060 0.2702
## 1500 meter buffer 0.0311900 -0.12760 0.2596
## 2250 meter buffer 0.0261500 -0.08821 0.2749
## 1750 meter buffer 0.0253800 -0.11130 0.2522
## 2000 meter buffer 0.0109100 -0.09835 0.2498
## 3000 meter buffer -0.0344400 -0.05023 0.2206
## 2500 meter buffer -0.0003667 -0.07148 0.2642
## 2750 meter buffer -0.0308900 -0.05897 0.2323
## 4000 meter buffer 0.0469100 -0.06100 0.2271
## 1000 meter buffer 0.0301900 -0.13170 0.2596
## 3250 meter buffer -0.0137100 -0.04806 0.2131
## 3500 meter buffer 0.0135300 -0.05270 0.2203
## 3750 meter buffer 0.0342300 -0.05710 0.2199
## 750 meter buffer 0.0875600 -0.13410 0.2508
## 500 meter buffer 0.0656500 -0.12500 0.1928
## 250 meter buffer 0.1312000 -0.02726 0.1534
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 5000 meter buffer -1.666e-02 0.009498 12 -460.069 945.6 0.00 0.444
## 4750 meter buffer -2.445e-02 0.016280 12 -460.272 946.0 0.41 0.363
## 4500 meter buffer -2.322e-02 0.036850 12 -461.227 947.9 2.32 0.139
## 4250 meter buffer -1.812e-02 0.045610 12 -463.712 952.8 7.29 0.012
## 1250 meter buffer 1.330e-02 0.096600 12 -463.717 952.9 7.30 0.012
## 1500 meter buffer 3.691e-02 0.098590 12 -464.041 953.5 7.95 0.008
## 2250 meter buffer 6.001e-02 0.096140 12 -464.439 954.3 8.74 0.006
## 1750 meter buffer 4.822e-02 0.094730 12 -464.715 954.9 9.29 0.004
## 2000 meter buffer 4.571e-02 0.095260 12 -465.172 955.8 10.21 0.003
## 3000 meter buffer 2.613e-02 0.096310 12 -465.491 956.4 10.84 0.002
## 2500 meter buffer 4.630e-02 0.094370 12 -465.546 956.5 10.95 0.002
## 2750 meter buffer 2.768e-02 0.090510 12 -465.949 957.3 11.76 0.001
## 4000 meter buffer -2.052e-03 0.059430 12 -466.033 957.5 11.93 0.001
## 1000 meter buffer -6.607e-05 0.137300 12 -466.066 957.6 11.99 0.001
## 3250 meter buffer 1.475e-02 0.080040 12 -466.272 958.0 12.41 0.001
## 3500 meter buffer 1.939e-02 0.081380 12 -466.715 958.9 13.29 0.001
## 3750 meter buffer 1.475e-02 0.078790 12 -466.789 959.0 13.44 0.001
## 750 meter buffer -3.534e-02 0.191500 12 -469.320 964.1 18.50 0.000
## 500 meter buffer -7.750e-02 0.204300 12 -485.654 996.7 51.17 0.000
## 250 meter buffer 3.498e-02 0.030640 12 -503.904 1033.2 87.67 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Note about re-run global models: when we re-ran the global models with dropped correlated variables our top model for coyote did change (2250-5000). As with caribou this result is also more similar to the top anthro model (4750m)
Model summary 5000m
Let’s get the model summary for this model
summary(coyote_mods$`5000 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(coyote, absent_coyote) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## scale(lc_grassland) + scale(lc_developed) + scale(lc_forest) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 944.1 985.5 -460.1 920.1 220
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.1189 0.3449
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.387405 0.150996 -9.188 < 2e-16 ***
## scale(harvest) -0.084684 0.094828 -0.893 0.371845
## scale(seismic_lines) 0.266850 0.077153 3.459 0.000543 ***
## scale(seismic_lines_3D) -0.071213 0.074950 -0.950 0.342037
## scale(trails) -0.016663 0.072397 -0.230 0.817965
## scale(wells) 0.009498 0.111273 0.085 0.931979
## scale(osm_industrial) 0.353667 0.064260 5.504 3.72e-08 ***
## scale(pipeline_transmission_lines) 0.113690 0.106001 1.073 0.283479
## scale(lc_grassland) -0.201566 0.091179 -2.211 0.027059 *
## scale(lc_developed) 0.372383 0.090102 4.133 3.58e-05 ***
## scale(lc_forest) -0.159951 0.070220 -2.278 0.022735 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Everything looks good here
Subset models
Add later if deemed necessary.
Fisher
Global models
fisher_mods <- osm_final_df_2021_2022 %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(fisher, absent_fisher) ~
# HFI
scale(harvest) +
# scale(pipeline) +
# scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
# scale(transmission_lines) +
# scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# VEG covariates in numerical order
scale(lc_grassland) +
# scale(lc_coniferous) +
# scale(lc_broadleaf) +
# scale(lc_mixed) +
scale(lc_developed) +
# scale(lc_shrub) +
scale(lc_forest) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(fisher_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_dvl))
## 1000 meter buffer -3.100 + -0.2428 0.5727
## 750 meter buffer -3.079 + -0.2584 0.5548
## 500 meter buffer -3.063 + -0.3150 0.5010
## 2500 meter buffer -3.014 + -0.1366 0.5125
## 1250 meter buffer -3.061 + -0.2543 0.5114
## 2750 meter buffer -3.006 + -0.1379 0.4877
## 2250 meter buffer -3.026 + -0.1485 0.5460
## 3000 meter buffer -2.994 + -0.1425 0.4720
## 3250 meter buffer -2.987 + -0.1566 0.4501
## 2000 meter buffer -3.029 + -0.1700 0.5262
## 3500 meter buffer -2.980 + -0.1802 0.4492
## 1500 meter buffer -3.032 + -0.2348 0.4638
## 4250 meter buffer -2.979 + -0.2076 0.4306
## 1750 meter buffer -3.024 + -0.2122 0.4537
## 4500 meter buffer -2.984 + -0.2266 0.4129
## 4000 meter buffer -2.974 + -0.1987 0.4324
## 3750 meter buffer -2.973 + -0.1958 0.4323
## 4750 meter buffer -2.990 + -0.2586 0.3990
## 5000 meter buffer -2.989 + -0.2811 0.3666
## 250 meter buffer -3.002 + -0.3663 0.2780
## cnd(scl(lc_frs)) cnd(scl(lc_grs)) cnd(scl(osm_ind))
## 1000 meter buffer 0.7062 0.21160 -0.020860
## 750 meter buffer 0.7434 0.22500 -0.057030
## 500 meter buffer 0.6974 0.26370 -0.195300
## 2500 meter buffer 0.5431 -0.11200 0.031300
## 1250 meter buffer 0.6093 0.15920 -0.044020
## 2750 meter buffer 0.5312 -0.16380 0.046910
## 2250 meter buffer 0.5793 -0.01139 0.037880
## 3000 meter buffer 0.4993 -0.21700 0.053230
## 3250 meter buffer 0.4830 -0.24140 0.001689
## 2000 meter buffer 0.5935 0.01941 0.027620
## 3500 meter buffer 0.4816 -0.21500 -0.055970
## 1500 meter buffer 0.5691 0.10920 -0.017430
## 4250 meter buffer 0.4670 -0.18050 -0.218500
## 1750 meter buffer 0.5801 0.05582 0.029350
## 4500 meter buffer 0.4764 -0.15340 -0.266600
## 4000 meter buffer 0.4705 -0.16090 -0.139500
## 3750 meter buffer 0.4690 -0.17660 -0.102500
## 4750 meter buffer 0.4828 -0.12930 -0.325400
## 5000 meter buffer 0.4737 -0.10190 -0.327200
## 250 meter buffer 0.4161 0.25530 -0.243500
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 1000 meter buffer -0.5916 -0.3746 2.617e-02
## 750 meter buffer -0.4845 -0.3410 1.612e-02
## 500 meter buffer -0.3401 -0.3199 -5.499e-02
## 2500 meter buffer -0.4934 -0.1681 -8.696e-02
## 1250 meter buffer -0.4646 -0.3741 -9.888e-03
## 2750 meter buffer -0.4607 -0.1506 -1.201e-01
## 2250 meter buffer -0.5188 -0.1975 -4.349e-02
## 3000 meter buffer -0.3993 -0.1326 -1.209e-01
## 3250 meter buffer -0.3517 -0.1280 -1.118e-01
## 2000 meter buffer -0.4857 -0.2378 -5.515e-02
## 3500 meter buffer -0.3216 -0.1315 -8.439e-02
## 1500 meter buffer -0.3868 -0.3431 -3.310e-02
## 4250 meter buffer -0.3576 -0.1193 -1.915e-02
## 1750 meter buffer -0.3844 -0.3047 -5.883e-02
## 4500 meter buffer -0.3674 -0.1200 -7.802e-03
## 4000 meter buffer -0.3664 -0.1227 -2.482e-02
## 3750 meter buffer -0.3284 -0.1345 -5.582e-02
## 4750 meter buffer -0.3377 -0.1396 6.660e-03
## 5000 meter buffer -0.2854 -0.1713 -9.999e-05
## 250 meter buffer -0.1913 -0.1936 -7.433e-02
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 1000 meter buffer 0.094360 0.04744 12 -269.750 564.9 0.00 0.526
## 750 meter buffer 0.065720 0.06888 12 -270.627 566.7 1.75 0.219
## 500 meter buffer 0.060860 0.09495 12 -271.319 568.1 3.14 0.109
## 2500 meter buffer -0.017210 0.18730 12 -272.236 569.9 4.97 0.044
## 1250 meter buffer 0.065770 0.04406 12 -272.682 570.8 5.86 0.028
## 2750 meter buffer -0.035120 0.21780 12 -272.778 571.0 6.06 0.025
## 2250 meter buffer -0.025000 0.12020 12 -273.366 572.2 7.23 0.014
## 3000 meter buffer -0.046220 0.22410 12 -273.682 572.8 7.86 0.010
## 3250 meter buffer -0.051490 0.26660 12 -273.895 573.2 8.29 0.008
## 2000 meter buffer -0.022470 0.10880 12 -273.934 573.3 8.37 0.008
## 3500 meter buffer -0.005697 0.25000 12 -275.380 576.2 11.26 0.002
## 1500 meter buffer 0.067490 0.03112 12 -275.814 577.1 12.13 0.001
## 4250 meter buffer 0.009936 0.34950 12 -275.904 577.2 12.31 0.001
## 1750 meter buffer 0.059910 0.06251 12 -276.090 577.6 12.68 0.001
## 4500 meter buffer 0.028680 0.37660 12 -276.236 577.9 12.97 0.001
## 4000 meter buffer 0.015420 0.30770 12 -276.305 578.0 13.11 0.001
## 3750 meter buffer 0.018250 0.27660 12 -276.378 578.2 13.26 0.001
## 4750 meter buffer 0.067160 0.37590 12 -276.558 578.5 13.62 0.001
## 5000 meter buffer 0.085460 0.33600 12 -277.564 580.6 15.63 0.000
## 250 meter buffer 0.008109 -0.09731 12 -278.248 581.9 17.00 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Note about re-run global models: when we re-ran the global models with dropped correlated variables our top model for fisher did not change but the weight did increase (0.266 - 0.526)
Model summary 1000m
Let’s print the summary for this model
summary(fisher_mods$`1000 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(fisher, absent_fisher) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## scale(lc_grassland) + scale(lc_developed) + scale(lc_forest) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 563.5 604.9 -269.7 539.5 220
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.2993 0.5471
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.09976 0.25568 -12.124 < 2e-16 ***
## scale(harvest) -0.24280 0.10744 -2.260 0.02383 *
## scale(seismic_lines) 0.02617 0.09586 0.273 0.78482
## scale(seismic_lines_3D) -0.37457 0.16863 -2.221 0.02633 *
## scale(trails) 0.09436 0.09191 1.027 0.30457
## scale(wells) 0.04744 0.11174 0.425 0.67117
## scale(osm_industrial) -0.02086 0.13613 -0.153 0.87821
## scale(pipeline_transmission_lines) -0.59158 0.19212 -3.079 0.00208 **
## scale(lc_grassland) 0.21162 0.11357 1.863 0.06240 .
## scale(lc_developed) 0.57274 0.12387 4.624 3.77e-06 ***
## scale(lc_forest) 0.70616 0.17333 4.074 4.62e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Everything looks good
Subset models
Add later if deemed necessary.
Grey wolf
Global models
wolf_mods <- osm_final_df_2021_2022 %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(grey_wolf, absent_grey_wolf) ~
# HFI
scale(harvest) +
# scale(pipeline) +
# scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
# scale(transmission_lines) +
# scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# VEG covariates in numerical order
scale(lc_grassland) +
# scale(lc_coniferous) +
# scale(lc_broadleaf) +
# scale(lc_mixed) +
scale(lc_developed) +
# scale(lc_shrub) +
scale(lc_forest) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(wolf_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_dvl))
## 3500 meter buffer -3.241 + 0.194200 0.0266400
## 3750 meter buffer -3.241 + 0.204400 0.0335100
## 3250 meter buffer -3.238 + 0.180700 0.0108700
## 4000 meter buffer -3.229 + 0.184100 -0.0353700
## 3000 meter buffer -3.238 + 0.188600 0.0003671
## 4250 meter buffer -3.210 + 0.122200 -0.0919800
## 2250 meter buffer -3.271 + 0.136400 0.3242000
## 2750 meter buffer -3.245 + 0.190300 0.0534600
## 2000 meter buffer -3.265 + 0.104400 0.3375000
## 2500 meter buffer -3.253 + 0.165100 0.1837000
## 4500 meter buffer -3.187 + -0.023640 -0.2649000
## 1750 meter buffer -3.248 + 0.077000 0.3189000
## 4750 meter buffer -3.191 + 0.076090 -0.0583900
## 5000 meter buffer -3.188 + 0.077810 -0.0061590
## 1500 meter buffer -3.224 + 0.072340 0.2963000
## 1250 meter buffer -3.190 + 0.081340 0.1546000
## 1000 meter buffer -3.173 + 0.067840 0.0476200
## 250 meter buffer -3.212 + 0.024560 -0.1363000
## 500 meter buffer -3.158 + -0.002271 -0.0650100
## 750 meter buffer -3.150 + 0.069940 -0.0863700
## cnd(scl(lc_frs)) cnd(scl(lc_grs)) cnd(scl(osm_ind))
## 3500 meter buffer 0.130800 0.56010 -0.18350
## 3750 meter buffer 0.140200 0.59060 -0.17140
## 3250 meter buffer 0.127500 0.49870 -0.17060
## 4000 meter buffer 0.145500 0.58870 -0.13710
## 3000 meter buffer 0.103100 0.43690 -0.17490
## 4250 meter buffer 0.144100 0.54810 -0.09837
## 2250 meter buffer 0.197800 0.26460 -0.33710
## 2750 meter buffer 0.120000 0.38450 -0.18850
## 2000 meter buffer 0.218100 0.21950 -0.37560
## 2500 meter buffer 0.141900 0.30710 -0.25880
## 4500 meter buffer 0.119900 0.47800 -0.04723
## 1750 meter buffer 0.206800 0.15510 -0.34620
## 4750 meter buffer 0.153800 0.51130 -0.08593
## 5000 meter buffer 0.155400 0.50950 -0.09069
## 1500 meter buffer 0.173800 0.09912 -0.32430
## 1250 meter buffer 0.093510 0.08397 -0.27580
## 1000 meter buffer 0.044090 0.06164 -0.19450
## 250 meter buffer -0.016460 -0.06205 0.07144
## 500 meter buffer -0.110400 -0.12390 -0.06536
## 750 meter buffer -0.008995 0.04169 -0.10910
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 3500 meter buffer -0.163300 -0.7269 -0.02457
## 3750 meter buffer -0.204500 -0.6843 -0.01182
## 3250 meter buffer -0.111800 -0.7746 -0.01252
## 4000 meter buffer -0.239500 -0.6796 -0.02869
## 3000 meter buffer -0.073020 -0.8350 0.00878
## 4250 meter buffer -0.237500 -0.7003 -0.03495
## 2250 meter buffer -0.016120 -0.8264 0.03814
## 2750 meter buffer -0.053470 -0.8663 0.02660
## 2000 meter buffer 0.002686 -0.8083 0.01924
## 2500 meter buffer -0.017180 -0.8666 0.04403
## 4500 meter buffer -0.162500 -0.8086 -0.05614
## 1750 meter buffer 0.007225 -0.7936 -0.06795
## 4750 meter buffer -0.306600 -0.6688 -0.04269
## 5000 meter buffer -0.359800 -0.6356 -0.03932
## 1500 meter buffer -0.000869 -0.7547 -0.04433
## 1250 meter buffer -0.006261 -0.7537 -0.00261
## 1000 meter buffer 0.008651 -0.7957 0.04657
## 250 meter buffer 0.028610 -1.1090 0.01309
## 500 meter buffer 0.016400 -0.8089 0.06742
## 750 meter buffer 0.037890 -0.8334 0.11680
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 3500 meter buffer 0.09244 -0.40670 12 -241.030 507.5 0.00 0.303
## 3750 meter buffer 0.07760 -0.43920 12 -241.253 507.9 0.45 0.242
## 3250 meter buffer 0.10370 -0.39190 12 -241.810 509.0 1.56 0.139
## 4000 meter buffer 0.06075 -0.37870 12 -241.976 509.4 1.89 0.118
## 3000 meter buffer 0.09972 -0.34470 12 -242.930 511.3 3.80 0.045
## 4250 meter buffer 0.07843 -0.33630 12 -243.276 512.0 4.49 0.032
## 2250 meter buffer 0.19990 -0.49620 12 -243.431 512.3 4.80 0.027
## 2750 meter buffer 0.11660 -0.33950 12 -243.538 512.5 5.02 0.025
## 2000 meter buffer 0.22960 -0.45730 12 -243.663 512.8 5.27 0.022
## 2500 meter buffer 0.16640 -0.40500 12 -243.982 513.4 5.90 0.016
## 4500 meter buffer 0.11890 -0.23370 12 -244.112 513.6 6.16 0.014
## 1750 meter buffer 0.22580 -0.42120 12 -244.635 514.7 7.21 0.008
## 4750 meter buffer 0.08722 -0.29140 12 -245.213 515.9 8.36 0.005
## 5000 meter buffer 0.08120 -0.29240 12 -245.842 517.1 9.62 0.002
## 1500 meter buffer 0.21080 -0.36820 12 -245.924 517.3 9.79 0.002
## 1250 meter buffer 0.15810 -0.28710 12 -247.955 521.3 13.85 0.000
## 1000 meter buffer 0.10350 -0.27580 12 -249.381 524.2 16.70 0.000
## 250 meter buffer 0.22050 -0.06846 12 -249.716 524.9 17.37 0.000
## 500 meter buffer 0.10090 -0.30240 12 -249.844 525.1 17.63 0.000
## 750 meter buffer 0.03225 -0.12850 12 -251.432 528.3 20.80 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Note about re-run global models: when we re-ran the global models with dropped correlated variables our top model for wolf did change (500 - 3500), as with coyote and caribou this is more similar to top anthro modelthan it was before AND actually the same as the landscape model
Model summary 3500m
Let’s get the model summary for this model
summary(wolf_mods$`3500 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(grey_wolf, absent_grey_wolf) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## scale(lc_grassland) + scale(lc_developed) + scale(lc_forest) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 506.1 547.4 -241.0 482.1 220
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.1303 0.3609
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.24108 0.19902 -16.285 < 2e-16 ***
## scale(harvest) 0.19422 0.16181 1.200 0.230006
## scale(seismic_lines) -0.02457 0.13065 -0.188 0.850840
## scale(seismic_lines_3D) -0.72687 0.29274 -2.483 0.013026 *
## scale(trails) 0.09244 0.11377 0.812 0.416513
## scale(wells) -0.40674 0.27482 -1.480 0.138866
## scale(osm_industrial) -0.18346 0.16334 -1.123 0.261349
## scale(pipeline_transmission_lines) -0.16326 0.19750 -0.827 0.408440
## scale(lc_grassland) 0.56008 0.15213 3.682 0.000232 ***
## scale(lc_developed) 0.02664 0.28028 0.095 0.924282
## scale(lc_forest) 0.13081 0.19835 0.659 0.509581
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Subset models
Add later if deemed necessary.
Lynx
Global models
lynx_mods <- osm_final_df_2021_2022 %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(lynx, absent_lynx) ~
# HFI
scale(harvest) +
# scale(pipeline) +
# scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
# scale(transmission_lines) +
# scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# VEG covariates in numerical order
scale(lc_grassland) +
# scale(lc_coniferous) +
# scale(lc_broadleaf) +
# scale(lc_mixed) +
scale(lc_developed) +
# scale(lc_shrub) +
scale(lc_forest) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(lynx_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_dvl))
## 250 meter buffer -1.898 + -0.302900 -0.13490
## 500 meter buffer -1.884 + -0.217500 -0.09147
## 3250 meter buffer -1.888 + 0.169200 0.14350
## 3000 meter buffer -1.887 + 0.159700 0.16200
## 1000 meter buffer -1.873 + -0.106800 -0.01345
## 3500 meter buffer -1.886 + 0.172100 0.13140
## 1250 meter buffer -1.873 + -0.089820 0.07802
## 750 meter buffer -1.875 + -0.150000 -0.04747
## 3750 meter buffer -1.885 + 0.175600 0.12880
## 2750 meter buffer -1.883 + 0.138700 0.16530
## 1500 meter buffer -1.873 + -0.059070 0.14370
## 4250 meter buffer -1.882 + 0.180100 0.09102
## 4000 meter buffer -1.883 + 0.177300 0.10850
## 5000 meter buffer -1.875 + 0.129400 0.04350
## 4750 meter buffer -1.875 + 0.145900 0.05866
## 4500 meter buffer -1.878 + 0.170800 0.06515
## 2000 meter buffer -1.873 + 0.040980 0.18330
## 2500 meter buffer -1.877 + 0.105100 0.18000
## 1750 meter buffer -1.870 + -0.006825 0.15910
## 2250 meter buffer -1.874 + 0.079930 0.18790
## cnd(scl(lc_frs)) cnd(scl(lc_grs)) cnd(scl(osm_ind))
## 250 meter buffer -0.0185600 0.01085 0.1005000
## 500 meter buffer 0.0469000 0.16350 0.1223000
## 3250 meter buffer 0.0121600 0.16580 0.0087360
## 3000 meter buffer 0.0203100 0.17250 0.0027110
## 1000 meter buffer 0.0776200 0.13000 0.0816700
## 3500 meter buffer -0.0003699 0.13620 0.0015910
## 1250 meter buffer 0.0878800 0.14980 0.0438800
## 750 meter buffer 0.0754100 0.12480 0.1076000
## 3750 meter buffer -0.0065780 0.09998 0.0004248
## 2750 meter buffer 0.0220000 0.16320 0.0039440
## 1500 meter buffer 0.0902000 0.14670 0.0195500
## 4250 meter buffer 0.0024930 0.06487 0.0250700
## 4000 meter buffer -0.0039010 0.07578 0.0103100
## 5000 meter buffer 0.0295100 -0.01200 0.0781900
## 4750 meter buffer 0.0196100 -0.01061 0.0629900
## 4500 meter buffer 0.0124200 0.02205 0.0472300
## 2000 meter buffer 0.0804900 0.10730 0.0145200
## 2500 meter buffer 0.0380200 0.11930 0.0075050
## 1750 meter buffer 0.0837900 0.10560 0.0182500
## 2250 meter buffer 0.0602100 0.09286 0.0140300
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 250 meter buffer 0.00373 0.07554 -0.172500
## 500 meter buffer -0.09004 0.07648 -0.090880
## 3250 meter buffer -0.10370 0.11650 -0.137800
## 3000 meter buffer -0.10520 0.12050 -0.106900
## 1000 meter buffer -0.03502 0.07623 -0.001316
## 3500 meter buffer -0.08795 0.10400 -0.158800
## 1250 meter buffer -0.05108 0.09034 0.061220
## 750 meter buffer -0.04502 0.07206 -0.038620
## 3750 meter buffer -0.07978 0.09587 -0.175600
## 2750 meter buffer -0.09281 0.11960 -0.083160
## 1500 meter buffer -0.05109 0.10160 0.067110
## 4250 meter buffer -0.07660 0.09533 -0.177500
## 4000 meter buffer -0.08150 0.09409 -0.175500
## 5000 meter buffer -0.07420 0.10170 -0.167300
## 4750 meter buffer -0.06444 0.09295 -0.163800
## 4500 meter buffer -0.07609 0.09271 -0.169700
## 2000 meter buffer -0.04631 0.10570 0.026900
## 2500 meter buffer -0.08202 0.11110 -0.062450
## 1750 meter buffer -0.03841 0.10150 0.043040
## 2250 meter buffer -0.06802 0.10680 -0.035860
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 250 meter buffer 0.04983 -0.2173 12 -436.043 897.5 0.00 0.982
## 500 meter buffer -0.01788 -0.2137 12 -440.225 905.9 8.36 0.015
## 3250 meter buffer 0.09311 -0.2843 12 -443.642 912.7 15.20 0.000
## 3000 meter buffer 0.09253 -0.3065 12 -443.650 912.7 15.21 0.000
## 1000 meter buffer 0.10100 -0.2316 12 -444.107 913.6 16.13 0.000
## 3500 meter buffer 0.09992 -0.2578 12 -444.112 913.6 16.14 0.000
## 1250 meter buffer 0.11410 -0.2692 12 -444.222 913.9 16.36 0.000
## 750 meter buffer 0.03278 -0.1957 12 -444.487 914.4 16.89 0.000
## 3750 meter buffer 0.09758 -0.2372 12 -444.527 914.5 16.97 0.000
## 2750 meter buffer 0.08370 -0.3161 12 -444.687 914.8 17.29 0.000
## 1500 meter buffer 0.12620 -0.2865 12 -444.871 915.2 17.66 0.000
## 4250 meter buffer 0.09542 -0.1957 12 -445.179 915.8 18.27 0.000
## 4000 meter buffer 0.09328 -0.2017 12 -445.338 916.1 18.59 0.000
## 5000 meter buffer 0.09053 -0.1283 12 -445.736 916.9 19.39 0.000
## 4750 meter buffer 0.09617 -0.1402 12 -445.792 917.0 19.50 0.000
## 4500 meter buffer 0.09432 -0.1481 12 -445.824 917.1 19.56 0.000
## 2000 meter buffer 0.11280 -0.3164 12 -445.995 917.4 19.90 0.000
## 2500 meter buffer 0.08585 -0.3024 12 -446.278 918.0 20.47 0.000
## 1750 meter buffer 0.11680 -0.2939 12 -446.313 918.1 20.54 0.000
## 2250 meter buffer 0.09530 -0.2986 12 -446.614 918.7 21.14 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Note about re-run global models: when we re-ran the global models with dropped correlated variables our top model for lynx did change (1250 - 250), as with others this is more similar to top anthro model (acutally the same as anthro model).
Model summary 250m
Let’s get the model summary
summary(lynx_mods$`250 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(lynx, absent_lynx) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## scale(lc_grassland) + scale(lc_developed) + scale(lc_forest) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 896.1 937.4 -436.0 872.1 220
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.1573 0.3967
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.89843 0.17291 -10.980 < 2e-16 ***
## scale(harvest) -0.30288 0.09009 -3.362 0.000774 ***
## scale(seismic_lines) -0.17254 0.07219 -2.390 0.016840 *
## scale(seismic_lines_3D) 0.07554 0.05476 1.379 0.167755
## scale(trails) 0.04983 0.05307 0.939 0.347748
## scale(wells) -0.21728 0.07372 -2.947 0.003204 **
## scale(osm_industrial) 0.10052 0.05392 1.864 0.062295 .
## scale(pipeline_transmission_lines) 0.00373 0.06477 0.058 0.954079
## scale(lc_grassland) 0.01085 0.07148 0.152 0.879349
## scale(lc_developed) -0.13488 0.07953 -1.696 0.089878 .
## scale(lc_forest) -0.01856 0.07703 -0.241 0.809563
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Everything looks good so far
Subset models
Add later if deemed necessary.
Moose
Global model
moose_mods <- osm_final_df_2021_2022 %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(moose, absent_moose) ~
# HFI
scale(harvest) +
# scale(pipeline) +
# scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
# scale(transmission_lines) +
# scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# VEG covariates in numerical order
scale(lc_grassland) +
# scale(lc_coniferous) +
# scale(lc_broadleaf) +
# scale(lc_mixed) +
scale(lc_developed) +
# scale(lc_shrub) +
scale(lc_forest) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(moose_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_dvl))
## 250 meter buffer -1.850 + -0.010620 -0.400100
## 500 meter buffer -1.845 + 0.007943 -0.247700
## 750 meter buffer -1.843 + 0.022290 -0.171400
## 1000 meter buffer -1.839 + 0.002889 -0.171700
## 1250 meter buffer -1.832 + 0.022190 -0.153200
## 1500 meter buffer -1.832 + 0.017690 -0.078240
## 1750 meter buffer -1.829 + 0.005968 -0.067750
## 2500 meter buffer -1.827 + -0.020630 -0.011380
## 2750 meter buffer -1.828 + -0.020800 0.003767
## 3000 meter buffer -1.827 + -0.024110 0.032780
## 3250 meter buffer -1.825 + -0.033290 0.047790
## 2000 meter buffer -1.827 + -0.010090 -0.018670
## 2250 meter buffer -1.826 + -0.008670 -0.025810
## 3500 meter buffer -1.823 + -0.046410 0.047440
## 3750 meter buffer -1.820 + -0.058910 0.046930
## 4000 meter buffer -1.818 + -0.067570 0.043240
## 5000 meter buffer -1.813 + -0.107100 0.047950
## 4250 meter buffer -1.815 + -0.080320 0.041380
## 4750 meter buffer -1.812 + -0.100500 0.062710
## 4500 meter buffer -1.813 + -0.092630 0.043460
## cnd(scl(lc_frs)) cnd(scl(lc_grs)) cnd(scl(osm_ind))
## 250 meter buffer -0.28910 -0.113100 -0.1808
## 500 meter buffer -0.20870 -0.028000 -0.2958
## 750 meter buffer -0.12680 0.012860 -0.3382
## 1000 meter buffer -0.15430 -0.005878 -0.3007
## 1250 meter buffer -0.16740 0.026720 -0.2424
## 1500 meter buffer -0.16390 0.026180 -0.2505
## 1750 meter buffer -0.15420 0.037200 -0.1993
## 2500 meter buffer -0.13540 -0.005974 -0.2374
## 2750 meter buffer -0.13330 0.007117 -0.2615
## 3000 meter buffer -0.12000 -0.015150 -0.2740
## 3250 meter buffer -0.10990 -0.029270 -0.2843
## 2000 meter buffer -0.13430 0.015520 -0.1879
## 2250 meter buffer -0.14130 -0.004633 -0.2027
## 3500 meter buffer -0.10790 -0.053200 -0.2895
## 3750 meter buffer -0.10970 -0.085060 -0.2964
## 4000 meter buffer -0.09742 -0.113900 -0.2822
## 5000 meter buffer -0.06937 -0.227500 -0.1979
## 4250 meter buffer -0.08018 -0.123200 -0.2578
## 4750 meter buffer -0.06349 -0.160100 -0.2124
## 4500 meter buffer -0.07208 -0.128100 -0.2348
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 250 meter buffer 0.04241 0.1464000 0.105100
## 500 meter buffer 0.07822 0.1164000 0.073920
## 750 meter buffer 0.10340 0.0916400 0.003156
## 1000 meter buffer 0.10280 0.0680300 -0.065830
## 1250 meter buffer 0.10130 0.0580600 -0.046040
## 1500 meter buffer 0.10060 0.0567600 -0.060180
## 1750 meter buffer 0.11300 0.0397600 -0.068110
## 2500 meter buffer 0.20240 -0.0171900 -0.014270
## 2750 meter buffer 0.21300 -0.0211100 -0.014690
## 3000 meter buffer 0.21860 -0.0239100 -0.032260
## 3250 meter buffer 0.21820 -0.0276000 -0.032760
## 2000 meter buffer 0.13050 0.0213800 -0.037160
## 2250 meter buffer 0.16860 -0.0006895 -0.022440
## 3500 meter buffer 0.22710 -0.0353600 -0.031970
## 3750 meter buffer 0.23510 -0.0442200 -0.028400
## 4000 meter buffer 0.22770 -0.0519300 -0.028000
## 5000 meter buffer 0.25780 -0.0746000 0.026540
## 4250 meter buffer 0.20750 -0.0498500 -0.014470
## 4750 meter buffer 0.21000 -0.0522300 0.020230
## 4500 meter buffer 0.19890 -0.0486300 0.003291
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 250 meter buffer 0.08051 -0.03034 12 -427.083 879.6 0.00 0.941
## 500 meter buffer 0.05705 -0.05982 12 -430.593 886.6 7.02 0.028
## 750 meter buffer 0.04281 -0.12910 12 -431.149 887.7 8.13 0.016
## 1000 meter buffer 0.06116 -0.16430 12 -431.351 888.1 8.54 0.013
## 1250 meter buffer 0.07610 -0.14410 12 -434.278 894.0 14.39 0.001
## 1500 meter buffer 0.08446 -0.19610 12 -434.971 895.4 15.78 0.000
## 1750 meter buffer 0.06579 -0.25990 12 -435.704 896.8 17.24 0.000
## 2500 meter buffer 0.08758 -0.27100 12 -436.444 898.3 18.72 0.000
## 2750 meter buffer 0.08737 -0.25670 12 -436.609 898.6 19.05 0.000
## 3000 meter buffer 0.09077 -0.26310 12 -436.675 898.8 19.18 0.000
## 3250 meter buffer 0.09238 -0.24350 12 -437.098 899.6 20.03 0.000
## 2000 meter buffer 0.06643 -0.29840 12 -437.156 899.7 20.14 0.000
## 2250 meter buffer 0.07148 -0.27640 12 -437.203 899.8 20.24 0.000
## 3500 meter buffer 0.09331 -0.21970 12 -437.285 900.0 20.40 0.000
## 3750 meter buffer 0.09740 -0.17820 12 -437.450 900.3 20.73 0.000
## 4000 meter buffer 0.10050 -0.15070 12 -437.796 901.0 21.42 0.000
## 5000 meter buffer 0.10800 -0.14820 12 -438.003 901.4 21.84 0.000
## 4250 meter buffer 0.10400 -0.13820 12 -438.659 902.7 23.15 0.000
## 4750 meter buffer 0.10470 -0.17150 12 -438.989 903.4 23.81 0.000
## 4500 meter buffer 0.10500 -0.14420 12 -439.210 903.8 24.25 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Note about re-run global models: when we re-ran the global models with dropped correlated variables our top model for moose did not change but the weight did increase (0.504 - 0.941). This top model is the same as the landscape model.
Model summary 250m
Let’s get the model summary
summary(moose_mods$`250 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(moose, absent_moose) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## scale(lc_grassland) + scale(lc_developed) + scale(lc_forest) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 878.2 919.5 -427.1 854.2 220
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.3014 0.549
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.84972 0.23284 -7.944 1.95e-15 ***
## scale(harvest) -0.01062 0.06092 -0.174 0.861546
## scale(seismic_lines) 0.10507 0.06305 1.666 0.095643 .
## scale(seismic_lines_3D) 0.14644 0.05514 2.656 0.007909 **
## scale(trails) 0.08051 0.05258 1.531 0.125766
## scale(wells) -0.03034 0.06728 -0.451 0.652026
## scale(osm_industrial) -0.18077 0.07146 -2.530 0.011420 *
## scale(pipeline_transmission_lines) 0.04241 0.06953 0.610 0.541879
## scale(lc_grassland) -0.11309 0.06932 -1.631 0.102817
## scale(lc_developed) -0.40009 0.08750 -4.572 4.82e-06 ***
## scale(lc_forest) -0.28907 0.08148 -3.548 0.000388 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Looks good for now
Subset models
Add later if deemed necessary
Red fox
Global models
fox_mods <- osm_final_df_2021_2022 %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(red_fox, absent_red_fox) ~
# HFI
scale(harvest) +
# scale(pipeline) +
# scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
# scale(transmission_lines) +
# scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# VEG covariates in numerical order
scale(lc_grassland) +
# scale(lc_coniferous) +
# scale(lc_broadleaf) +
# scale(lc_mixed) +
scale(lc_developed) +
# scale(lc_shrub) +
scale(lc_forest) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(fox_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_dvl))
## 4750 meter buffer -4.711 + 1.719e-01 0.3922
## 5000 meter buffer -4.438 + 1.962e-01 0.6080
## 4500 meter buffer -4.705 + 1.480e-01 0.3742
## 1500 meter buffer -4.533 + -2.459e-01 0.4051
## 1750 meter buffer -4.593 + -1.465e-01 0.2375
## 4250 meter buffer -4.582 + 1.112e-01 0.2584
## 4000 meter buffer -4.307 + 8.981e-02 0.5899
## 1250 meter buffer -4.413 + -2.390e-01 0.4211
## 3750 meter buffer -4.259 + 3.255e-02 0.6150
## 2000 meter buffer -4.463 + -9.591e-02 0.3378
## 2250 meter buffer -4.206 + -5.169e-02 0.7396
## 2500 meter buffer -4.185 + -2.673e-03 0.7367
## 1000 meter buffer -4.423 + -2.050e-01 0.5976
## 2750 meter buffer -4.166 + 4.747e-02 0.6695
## 500 meter buffer -4.332 + -1.192e-01 0.6214
## 3500 meter buffer -4.159 + 3.266e-05 0.6637
## 3000 meter buffer -4.150 + 5.339e-02 0.6575
## 3250 meter buffer -4.148 + 3.464e-02 0.5886
## 750 meter buffer -4.350 + -2.257e-01 0.4363
## 250 meter buffer -4.314 + -7.559e-02 0.5478
## cnd(scl(lc_frs)) cnd(scl(lc_grs)) cnd(scl(osm_ind))
## 4750 meter buffer -0.79370 0.4439 0.211500
## 5000 meter buffer -0.62830 0.4975 0.024190
## 4500 meter buffer -0.79680 0.3944 0.204900
## 1500 meter buffer -0.12630 0.5745 0.167300
## 1750 meter buffer -0.16780 0.5893 0.303300
## 4250 meter buffer -0.76450 0.2999 0.039530
## 4000 meter buffer -0.45660 0.4669 -0.497800
## 1250 meter buffer -0.10870 0.5177 0.165400
## 3750 meter buffer -0.31440 0.5524 -0.567400
## 2000 meter buffer -0.10970 0.6124 0.234700
## 2250 meter buffer -0.01560 0.7021 -0.006098
## 2500 meter buffer -0.06085 0.6474 -0.019590
## 1000 meter buffer -0.04114 0.5440 0.130100
## 2750 meter buffer -0.11750 0.5551 -0.053910
## 500 meter buffer 0.36550 0.3196 0.363800
## 3500 meter buffer -0.20010 0.5830 -0.253700
## 3000 meter buffer -0.13370 0.5417 -0.108200
## 3250 meter buffer -0.13510 0.5875 -0.139000
## 750 meter buffer 0.11420 0.5122 0.221500
## 250 meter buffer 0.14060 0.2998 0.192600
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 4750 meter buffer 0.85130 -0.64070 -0.566600
## 5000 meter buffer 0.36200 -1.29200 -0.236500
## 4500 meter buffer 0.83390 -0.59590 -0.517800
## 1500 meter buffer -0.26360 -0.38790 -0.121800
## 1750 meter buffer -0.09303 -0.47180 -0.191300
## 4250 meter buffer 0.75990 -0.72240 -0.414300
## 4000 meter buffer 0.20960 -1.01300 -0.054080
## 1250 meter buffer -0.58180 -0.39440 -0.092140
## 3750 meter buffer 0.10940 -0.89170 -0.001071
## 2000 meter buffer -0.22120 -0.50570 -0.141000
## 2250 meter buffer -0.47670 -0.65920 0.162900
## 2500 meter buffer -0.39530 -0.62460 0.136700
## 1000 meter buffer -0.92470 -0.27190 -0.067030
## 2750 meter buffer -0.29320 -0.65570 0.110900
## 500 meter buffer -0.34870 -0.25270 -0.173400
## 3500 meter buffer -0.05723 -0.69820 0.040430
## 3000 meter buffer -0.18210 -0.64930 0.102700
## 3250 meter buffer -0.17200 -0.70830 0.088240
## 750 meter buffer -0.37060 -0.16540 -0.052970
## 250 meter buffer -0.40470 0.04142 -0.087320
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 4750 meter buffer -0.8161 -0.8512 12 -163.481 352.4 0.00 0.367
## 5000 meter buffer -0.8335 -0.4373 12 -163.494 352.4 0.03 0.362
## 4500 meter buffer -0.8812 -0.8091 12 -164.121 353.7 1.28 0.193
## 1500 meter buffer -0.7754 -0.5035 12 -165.990 357.4 5.02 0.030
## 1750 meter buffer -0.8772 -0.4425 12 -166.518 358.5 6.08 0.018
## 4250 meter buffer -0.7284 -0.4446 12 -166.799 359.0 6.64 0.013
## 4000 meter buffer -0.6360 -0.1574 12 -167.584 360.6 8.21 0.006
## 1250 meter buffer -0.7080 -0.1627 12 -168.567 362.6 10.17 0.002
## 3750 meter buffer -0.4702 -0.1406 12 -168.618 362.7 10.28 0.002
## 2000 meter buffer -0.7710 -0.2948 12 -168.769 363.0 10.58 0.002
## 2250 meter buffer -0.6500 -0.1999 12 -168.871 363.2 10.78 0.002
## 2500 meter buffer -0.6585 -0.2297 12 -169.197 363.8 11.43 0.001
## 1000 meter buffer -0.6733 -0.1584 12 -169.244 363.9 11.53 0.001
## 2750 meter buffer -0.6457 -0.1508 12 -170.527 366.5 14.09 0.000
## 500 meter buffer -0.6186 0.2243 12 -170.872 367.2 14.78 0.000
## 3500 meter buffer -0.4202 -0.2561 12 -170.935 367.3 14.91 0.000
## 3000 meter buffer -0.5830 -0.1700 12 -171.405 368.2 15.85 0.000
## 3250 meter buffer -0.4497 -0.1065 12 -171.866 369.2 16.77 0.000
## 750 meter buffer -0.6205 -0.1641 12 -172.440 370.3 17.92 0.000
## 250 meter buffer -0.8329 0.1147 12 -173.426 372.3 19.89 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Note about re-run global models: when we re-ran the global models with dropped correlated variables our top model for fox did not change but the weight dropped (0.855 - 0.367). This top model is the same as the landscape model.
Model summary 4750m
Let’s get the model summary
summary(fox_mods$`4750 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(red_fox, absent_red_fox) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## scale(lc_grassland) + scale(lc_developed) + scale(lc_forest) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 351.0 392.3 -163.5 327.0 220
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 1.334 1.155
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.7107 0.6049 -7.788 6.83e-15 ***
## scale(harvest) 0.1719 0.2779 0.619 0.53603
## scale(seismic_lines) -0.5666 0.2969 -1.909 0.05632 .
## scale(seismic_lines_3D) -0.6407 0.6589 -0.972 0.33086
## scale(trails) -0.8161 0.4000 -2.040 0.04131 *
## scale(wells) -0.8512 0.4063 -2.095 0.03617 *
## scale(osm_industrial) 0.2115 0.3085 0.686 0.49298
## scale(pipeline_transmission_lines) 0.8513 0.4980 1.710 0.08736 .
## scale(lc_grassland) 0.4439 0.2471 1.796 0.07242 .
## scale(lc_developed) 0.3922 0.3527 1.112 0.26617
## scale(lc_forest) -0.7937 0.2942 -2.698 0.00697 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Looks good so far
Subset Models
Add later if deemed necessary.
White-tailed deer
Global models
deer_mods <- osm_final_df_2021_2022 %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
# have to include the `` around the white-tailed_deer or R won't recognize it as a variable because of the -
glmmTMB::glmmTMB(cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~
# HFI
scale(harvest) +
# scale(pipeline) +
# scale(roads) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
# scale(transmission_lines) +
# scale(veg_edges) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# VEG covariates in numerical order
scale(lc_grassland) +
# scale(lc_coniferous) +
# scale(lc_broadleaf) +
# scale(lc_mixed) +
scale(lc_developed) +
# scale(lc_shrub) +
scale(lc_forest) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(deer_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_dvl))
## 1500 meter buffer -0.2118 + -0.08698 -0.183400
## 4000 meter buffer -0.2228 + -0.20500 -0.087810
## 4500 meter buffer -0.2265 + -0.23140 -0.031440
## 4250 meter buffer -0.2245 + -0.22030 -0.066400
## 3750 meter buffer -0.2192 + -0.17860 -0.110400
## 4750 meter buffer -0.2282 + -0.25380 0.022210
## 3500 meter buffer -0.2173 + -0.15910 -0.117100
## 5000 meter buffer -0.2305 + -0.27600 0.043030
## 1250 meter buffer -0.2057 + -0.08675 -0.126400
## 1750 meter buffer -0.2068 + -0.09977 -0.172600
## 3250 meter buffer -0.2117 + -0.13330 -0.129100
## 3000 meter buffer -0.2055 + -0.11850 -0.120900
## 2000 meter buffer -0.2058 + -0.10150 -0.151500
## 2250 meter buffer -0.2033 + -0.11370 -0.137400
## 2750 meter buffer -0.2033 + -0.10730 -0.120900
## 2500 meter buffer -0.2043 + -0.10770 -0.123100
## 1000 meter buffer -0.1896 + -0.11220 -0.016110
## 750 meter buffer -0.1781 + -0.07972 0.033500
## 500 meter buffer -0.1829 + -0.01120 0.019780
## 250 meter buffer -0.1844 + 0.04015 -0.004757
## cnd(scl(lc_frs)) cnd(scl(lc_grs)) cnd(scl(osm_ind))
## 1500 meter buffer -0.06970 0.14370 0.3570
## 4000 meter buffer -0.08854 0.30650 0.2043
## 4500 meter buffer -0.09763 0.35500 0.1960
## 4250 meter buffer -0.09771 0.32570 0.1994
## 3750 meter buffer -0.09364 0.27870 0.2060
## 4750 meter buffer -0.08592 0.35770 0.1891
## 3500 meter buffer -0.09229 0.26600 0.2097
## 5000 meter buffer -0.08801 0.35140 0.1823
## 1250 meter buffer -0.05514 0.13260 0.3666
## 1750 meter buffer -0.07167 0.14090 0.3224
## 3250 meter buffer -0.09847 0.25220 0.2181
## 3000 meter buffer -0.10010 0.22910 0.2273
## 2000 meter buffer -0.07911 0.15440 0.2878
## 2250 meter buffer -0.08353 0.16050 0.2665
## 2750 meter buffer -0.09510 0.20520 0.2356
## 2500 meter buffer -0.09081 0.18920 0.2471
## 1000 meter buffer -0.01413 0.08257 0.3389
## 750 meter buffer -0.01587 0.05411 0.3300
## 500 meter buffer -0.05614 0.05070 0.3136
## 250 meter buffer -0.11240 0.02677 0.2171
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 1500 meter buffer 0.046930 -0.4547 -0.08996
## 4000 meter buffer -0.086430 -0.3458 -0.04042
## 4500 meter buffer -0.148500 -0.3000 -0.04713
## 4250 meter buffer -0.113200 -0.3221 -0.05741
## 3750 meter buffer -0.039670 -0.3672 -0.03410
## 4750 meter buffer -0.203300 -0.2724 -0.02822
## 3500 meter buffer -0.034270 -0.3659 -0.03981
## 5000 meter buffer -0.253600 -0.2489 -0.02063
## 1250 meter buffer 0.036410 -0.4348 -0.07624
## 1750 meter buffer 0.059220 -0.4230 -0.10160
## 3250 meter buffer -0.011200 -0.3568 -0.02783
## 3000 meter buffer -0.001411 -0.3422 -0.01609
## 2000 meter buffer 0.039230 -0.4015 -0.09493
## 2250 meter buffer 0.020150 -0.3672 -0.05275
## 2750 meter buffer 0.003639 -0.3422 -0.02911
## 2500 meter buffer 0.008484 -0.3450 -0.04237
## 1000 meter buffer 0.052860 -0.3421 -0.06818
## 750 meter buffer 0.098930 -0.2815 -0.01457
## 500 meter buffer 0.116700 -0.2645 0.03645
## 250 meter buffer 0.093550 -0.2163 -0.02893
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 1500 meter buffer 0.16850 0.35760 12 -559.536 1144.5 0.00 0.234
## 4000 meter buffer 0.22180 0.35280 12 -559.614 1144.7 0.16 0.216
## 4500 meter buffer 0.21000 0.31480 12 -559.740 1144.9 0.41 0.190
## 4250 meter buffer 0.21310 0.33550 12 -559.845 1145.1 0.62 0.172
## 3750 meter buffer 0.22320 0.36810 12 -560.452 1146.3 1.83 0.094
## 4750 meter buffer 0.20940 0.30110 12 -561.608 1148.6 4.15 0.029
## 3500 meter buffer 0.21630 0.37020 12 -562.094 1149.6 5.12 0.018
## 5000 meter buffer 0.21270 0.31850 12 -562.384 1150.2 5.70 0.014
## 1250 meter buffer 0.15510 0.31900 12 -562.447 1150.3 5.82 0.013
## 1750 meter buffer 0.16780 0.34860 12 -562.592 1150.6 6.11 0.011
## 3250 meter buffer 0.20690 0.37930 12 -562.987 1151.4 6.90 0.007
## 3000 meter buffer 0.18990 0.38040 12 -565.024 1155.5 10.98 0.001
## 2000 meter buffer 0.16500 0.35380 12 -565.365 1156.2 11.66 0.001
## 2250 meter buffer 0.17910 0.37890 12 -565.925 1157.3 12.78 0.000
## 2750 meter buffer 0.18800 0.38340 12 -566.039 1157.5 13.01 0.000
## 2500 meter buffer 0.18780 0.36990 12 -566.387 1158.2 13.70 0.000
## 1000 meter buffer 0.14890 0.27620 12 -572.128 1169.7 25.18 0.000
## 750 meter buffer 0.08956 0.25600 12 -583.589 1192.6 48.11 0.000
## 500 meter buffer 0.10640 0.24460 12 -588.883 1203.2 58.69 0.000
## 250 meter buffer 0.19490 0.06153 12 -595.626 1216.7 72.18 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Note about re-run global models: when we re-ran the global models with dropped correlated variables our top model for deer did not change but the weight dropped (0.974 - 0.234). This top model is the same as the landscape model.
Model summary 1500m
Let’s get the model summary
summary(deer_mods$`1500 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~ scale(harvest) +
## scale(seismic_lines) + scale(seismic_lines_3D) + scale(trails) +
## scale(wells) + scale(osm_industrial) + scale(pipeline_transmission_lines) +
## scale(lc_grassland) + scale(lc_developed) + scale(lc_forest) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 1143.1 1184.4 -559.5 1119.1 220
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 1.974 1.405
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.21180 0.57592 -0.368 0.71305
## scale(harvest) -0.08698 0.06553 -1.327 0.18440
## scale(seismic_lines) -0.08996 0.06637 -1.355 0.17530
## scale(seismic_lines_3D) -0.45472 0.09650 -4.712 2.45e-06 ***
## scale(trails) 0.16852 0.05189 3.247 0.00117 **
## scale(wells) 0.35761 0.07963 4.491 7.10e-06 ***
## scale(osm_industrial) 0.35696 0.05902 6.048 1.47e-09 ***
## scale(pipeline_transmission_lines) 0.04693 0.08430 0.557 0.57771
## scale(lc_grassland) 0.14373 0.06313 2.277 0.02281 *
## scale(lc_developed) -0.18344 0.08172 -2.245 0.02479 *
## scale(lc_forest) -0.06970 0.05879 -1.186 0.23578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Looks okay for now
Subset models
Add later if deemed necessary.
Top buffers
For convenience here is a list of the top buffer widths for each species
Initial run
* Black bear - 250m
* Caribou - 3000m
* Coyote - 2250m
* Fisher - 1000m
* Grey wolf - 500m
* Lynx - 1250m
* Moose - 250m
* Red fox - 4750
* White-tailed deer - 1500m
Re-run
* Black bear - 250m
* Caribou - 1000m
* Coyote - 5000m
* Fisher - 1000m
* Grey wolf - 3500m
* Lynx - 250m
* Moose - 250m
* Red fox - 4750m
* White-tailed deer - 1500m
Analysis Prep
Summary
While the global models are interesting and yield relevant results there is an issue of addressing autocorrelation with such large models that we don’t have the capacity to address at each spatial scale while still keeping the models the same so the only thing being compared is the buffer size.
It’s also relevant to wonder if these results would be different depending on if we look at exclusively anthropogenic or landscape features in a model and by reducing the variables in each model we can address this issue of autocorrelation much easier. So the next analysis will repeat what we’ve done above for each species but instead of using one global model for everything we will divide the data into disturbance and landscape covariates and run two separate models per species and see if the results differ.
Data prep
First let’s subset the data into the two new groups we will use for this analysis so we can check for auto correlations within the data sets
Anthropogenic
First let’s generate the developed data set which includes only the covariates for the human factors indices
osm_anthro_df_2021_2022 <- osm_final_df_2021_2022 %>%
# use purr to apply following data manipulation steps to all the buffer data frames
purrr::map(
~.x %>%
# de-select columns for landscape data (e.g. those that have the prefix 'lc_')
select(!starts_with('lc_'))
)Now let’s explore the correlations between variables at each buffer scale and see what features we may need to remove or group together (if appropriate)
Correlation plots
Now we need to make correlation plots for each buffer width to see
what variables are correlated at a given spatial scale. We can use
purrr::map() with the chart.Correlation()
function from the PerformanceAnalytics package to make
correlation plots with a specified method (e.g., pearson, spearman,
etc.) That also show histograms and scatterplots of each variable.
Let’s do this for the anthropogenic data first
osm_anthro_df_2021_2022 %>%
purrr::map(
~.x %>%
# select only columns with covaraites not other info
select(harvest:osm_industrial) %>%
# use chart.correlation
chart.Correlation(.,
histogram = TRUE,
method = "pearson")
)
## $`250 meter buffer`
## NULL
##
## $`500 meter buffer`
## NULL
##
## $`750 meter buffer`
## NULL
##
## $`1000 meter buffer`
## NULL
##
## $`1250 meter buffer`
## NULL
##
## $`1500 meter buffer`
## NULL
##
## $`1750 meter buffer`
## NULL
##
## $`2000 meter buffer`
## NULL
##
## $`2250 meter buffer`
## NULL
##
## $`2500 meter buffer`
## NULL
##
## $`2750 meter buffer`
## NULL
##
## $`3000 meter buffer`
## NULL
##
## $`3250 meter buffer`
## NULL
##
## $`3500 meter buffer`
## NULL
##
## $`3750 meter buffer`
## NULL
##
## $`4000 meter buffer`
## NULL
##
## $`4250 meter buffer`
## NULL
##
## $`4500 meter buffer`
## NULL
##
## $`4750 meter buffer`
## NULL
##
## $`5000 meter buffer`
## NULLThis creates a lot of plots (1 for each buffer size) and I haven’t found a way to label them by the buffer size yet unfortunately. But hopefully there will be trends among the buffer sizes to simplify the process of choosing covariates
in Rstudio you can click on the white square icon in the upper right-hand corner of the figures to open a new window with the plots so you can see the values easier
In order of buffer size I’m going to summarize covariates that are correlated below with their respective r2 values. I’m only listing correlations above 0.6 for each buffer size
250m roads + veg edges 0.71 this is to be expected as veg edges occur along the sides of roads (primarily) as well as other disturbance features. We will choose roads in the analysis over veg edges
500m pipeline + transmission lines 0.63
we may be able to combine these roads + veg edges 0.73750m
pipeline + transmission lines 0.70
roads + veg edges 0.741000m
pipeline + transmission lines 0.73
roads + veg edges 0.791250m
pipeline + transmission lines 0.72
roads + veg edges 0.811500m
pipeline + transmission lines 0.72
roads + veg edges 0.84
roads + wells 0.641750m pipeline + transmission lines 0.72
roads + veg edges 0.85
roads + wells 0.702000m
pipeline + transmission lines 0.72
roads + veg edges 0.86
roads + wells 0.732250m
pipeline + transmission lines 0.71
roads + veg edges 0.87
roads + wells 0.762500m
pipleine + roads 0.61
pipeline + transmission lines 0.70
roads + veg edges 0.88
roads + wells 0.76
roads + osm_industrial 0.612750m
pipeline + roads 0.61
pipeline + transmission lines 0.70
roads + veg edges 0.88
roads + wells 0.76
roads + osm_industrial 0.613000m
pipeline + roads 0.63
pipeline + transmission lines 0.69
roads + veg edges 0.88
roads + wells 0.77
roads + osm_industrial 0.63
at this point there are some obvious trends between roads and several features and pipelines and transmission lines. for simplicity sake I am only going to report correlations for variables not correlated with roads, as we will have to drop this covariate, and new pairs of highly correlated variables. I will also not specify the buffer size if it’s found in more than one buffer
- 3250m +
veg edges + wells 0.61
pipeline + wells 0.67
pipeline + wells 0.69
In summary we should merge pipelines and transmission lines, remove roads, and remove veg edges from the analysis to insure we don’t majorly violate assumptions of independence for our models. We will re-assess correlations after making these changes.
List of final variables
I first pasted the full list here and then deleted or added combined
variables as I went through each buffer’s correlation plot
harvest
seismic lines
seismic lines 3D
trails
pipeline and transmission lines
wells
osm_industrial
Reformat data
We’ll reformat the anthropogenicsubset data here and re-run the correlation plots
osm_anthro_df_2021_2022 <- osm_anthro_df_2021_2022 %>%
# use purrr
purrr::map(
~.x %>%
# mutate data to combine pipeline and transmission line
mutate(pipeline_transmission_lines = (pipeline + transmission_lines)) %>%
# remove correlated variables we won't include in subset models
select(! c(pipeline, roads, veg_edges, transmission_lines)) %>%
# relocate column so new covariate is with others
relocate(pipeline_transmission_lines,
.after = wells)
)Re-run correlation plots
Now we should double check that this solves any issues of highly correlated variables by re-running the plots from before
osm_anthro_df_2021_2022 %>%
purrr::map(
~.x %>%
# select only columns with covariates not other info
select(harvest:osm_industrial) %>%
# use chart.correlation
chart.Correlation(.,
histogram = TRUE,
method = "pearson")
)
## $`250 meter buffer`
## NULL
##
## $`500 meter buffer`
## NULL
##
## $`750 meter buffer`
## NULL
##
## $`1000 meter buffer`
## NULL
##
## $`1250 meter buffer`
## NULL
##
## $`1500 meter buffer`
## NULL
##
## $`1750 meter buffer`
## NULL
##
## $`2000 meter buffer`
## NULL
##
## $`2250 meter buffer`
## NULL
##
## $`2500 meter buffer`
## NULL
##
## $`2750 meter buffer`
## NULL
##
## $`3000 meter buffer`
## NULL
##
## $`3250 meter buffer`
## NULL
##
## $`3500 meter buffer`
## NULL
##
## $`3750 meter buffer`
## NULL
##
## $`4000 meter buffer`
## NULL
##
## $`4250 meter buffer`
## NULL
##
## $`4500 meter buffer`
## NULL
##
## $`4750 meter buffer`
## NULL
##
## $`5000 meter buffer`
## NULLWells and the new pipelines and transmission lines variable seem to become increasingly correlated as the buffer size increases but the max r squared value is 0.66 so we will leave it in for now and discuss with others about whether or not to remove it since it’s on the cusp.
Landscape
Now let’s do the same for the landscape data
osm_landscape_df_2021_2022 <- osm_final_df_2021_2022 %>%
# use purr to apply following data manipulation steps to all the buffer data frames
purrr::map(
~.x %>%
# de-select columns for developed data
select(!harvest:wells &
!osm_industrial)
)Now we will follow the same process with the landscape data for generating correlation plots and assessing whaich variables need to be removed or combined
Correlation plots
osm_landscape_df_2021_2022 %>%
purrr::map(
~.x %>%
# select only columns with covariates not other info
select(lc_grassland:lc_shrub) %>%
# use chart.correlation
chart.Correlation(.,
histogram = TRUE,
method = "pearson")
)
## $`250 meter buffer`
## NULL
##
## $`500 meter buffer`
## NULL
##
## $`750 meter buffer`
## NULL
##
## $`1000 meter buffer`
## NULL
##
## $`1250 meter buffer`
## NULL
##
## $`1500 meter buffer`
## NULL
##
## $`1750 meter buffer`
## NULL
##
## $`2000 meter buffer`
## NULL
##
## $`2250 meter buffer`
## NULL
##
## $`2500 meter buffer`
## NULL
##
## $`2750 meter buffer`
## NULL
##
## $`3000 meter buffer`
## NULL
##
## $`3250 meter buffer`
## NULL
##
## $`3500 meter buffer`
## NULL
##
## $`3750 meter buffer`
## NULL
##
## $`4000 meter buffer`
## NULL
##
## $`4250 meter buffer`
## NULL
##
## $`4500 meter buffer`
## NULL
##
## $`4750 meter buffer`
## NULL
##
## $`5000 meter buffer`
## NULLBelow is a summary of the different buffer sizes and correlated variables
250m
None >0.6500m coniferous + broadleaf 0.63
750m
coniferous + broadleaf 0.661000m
coniferous + braodleaf 0.68
Coniferous and broadleaf appear to be the only variables highly correlated so I think what we will do is combine all 3 forest types into one category, because we really aren’t interested in interpreting these individual variables but rather seeing if there are differences in which buffer is the best fit between development and landscape features.
Reformat data
osm_landscape_df_2021_2022 <- osm_landscape_df_2021_2022 %>%
# use purrr to apply the functions to all data frames
purrr::map(
~.x %>%
# use mutate to combine forest types into one variable
mutate(lc_forest = lc_coniferous + lc_broadleaf + lc_mixed) %>%
# remove old variables
select(! c(lc_coniferous, lc_broadleaf, lc_mixed)) %>%
# relocate the new column with. the other covariates
relocate(lc_forest,
.after = buff_dist)
)Looks good! Now let’s re-run the correlation plots
Re-run correlation plots
osm_landscape_df_2021_2022 %>%
purrr::map(
~.x %>%
# select only columns with covariates not other info
select(lc_forest:lc_shrub) %>%
# use chart.correlation
chart.Correlation(.,
histogram = TRUE,
method = "pearson")
)
## $`250 meter buffer`
## NULL
##
## $`500 meter buffer`
## NULL
##
## $`750 meter buffer`
## NULL
##
## $`1000 meter buffer`
## NULL
##
## $`1250 meter buffer`
## NULL
##
## $`1500 meter buffer`
## NULL
##
## $`1750 meter buffer`
## NULL
##
## $`2000 meter buffer`
## NULL
##
## $`2250 meter buffer`
## NULL
##
## $`2500 meter buffer`
## NULL
##
## $`2750 meter buffer`
## NULL
##
## $`3000 meter buffer`
## NULL
##
## $`3250 meter buffer`
## NULL
##
## $`3500 meter buffer`
## NULL
##
## $`3750 meter buffer`
## NULL
##
## $`4000 meter buffer`
## NULL
##
## $`4250 meter buffer`
## NULL
##
## $`4500 meter buffer`
## NULL
##
## $`4750 meter buffer`
## NULL
##
## $`5000 meter buffer`
## NULLWell crap. Turns out the primary landscape types are forest and shrub so when we combine all three forest types we now have an issue with autocorrelation with the shrub landcover type. Might need to rethink the approach here.
For now will use one or the other
Subset Analysis
Now that the two separate data files are ready we can run an analysis with each, similar to what we did with the global models and see if this changes our findings from the global model or if there are differences in results between the two model subsets
Anthropogenic
We will do another analysis for each species as we did previously.
Black bear
bear_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
# HFI that aren't correlated
scale(harvest) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))We will use the model.sel() function from the
MuMIn package to compare the global models for each buffer
width and see which buffer fits the black bear data best
model.sel(bear_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(osm_ind))
## 4000 meter buffer -0.6015 + 0.05542 -0.138100
## 3750 meter buffer -0.6014 + 0.05112 -0.139500
## 3500 meter buffer -0.6007 + 0.04791 -0.125600
## 4250 meter buffer -0.6007 + 0.06464 -0.132400
## 3250 meter buffer -0.6004 + 0.05222 -0.114300
## 250 meter buffer -0.5947 + 0.01448 0.040710
## 4500 meter buffer -0.6007 + 0.06538 -0.126200
## 4750 meter buffer -0.6005 + 0.06567 -0.121700
## 3000 meter buffer -0.5993 + 0.05158 -0.100700
## 2750 meter buffer -0.5987 + 0.05077 -0.087580
## 5000 meter buffer -0.5999 + 0.06688 -0.118400
## 2500 meter buffer -0.5977 + 0.04560 -0.074670
## 1750 meter buffer -0.5974 + 0.04413 -0.040290
## 2250 meter buffer -0.5971 + 0.04060 -0.063580
## 2000 meter buffer -0.5972 + 0.04329 -0.050380
## 1500 meter buffer -0.5971 + 0.04253 -0.027990
## 1000 meter buffer -0.5946 + 0.02595 -0.018970
## 1250 meter buffer -0.5956 + 0.03524 -0.008795
## 500 meter buffer -0.5941 + 0.03241 0.011470
## 750 meter buffer -0.5927 + 0.03608 0.004609
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 4000 meter buffer -0.16130 -0.08229 -0.12770
## 3750 meter buffer -0.15190 -0.08646 -0.11550
## 3500 meter buffer -0.15600 -0.08426 -0.10020
## 4250 meter buffer -0.15000 -0.08577 -0.14670
## 3250 meter buffer -0.16170 -0.08255 -0.09434
## 250 meter buffer -0.07483 -0.13580 -0.01959
## 4500 meter buffer -0.14820 -0.08715 -0.14890
## 4750 meter buffer -0.14590 -0.08946 -0.15040
## 3000 meter buffer -0.16370 -0.08422 -0.08098
## 2750 meter buffer -0.16570 -0.08585 -0.08105
## 5000 meter buffer -0.14450 -0.09169 -0.15030
## 2500 meter buffer -0.15460 -0.08323 -0.07427
## 1750 meter buffer -0.14610 -0.09229 -0.10730
## 2250 meter buffer -0.14930 -0.08317 -0.07225
## 2000 meter buffer -0.14760 -0.08610 -0.09301
## 1500 meter buffer -0.13530 -0.09867 -0.08145
## 1000 meter buffer -0.12710 -0.11110 -0.04637
## 1250 meter buffer -0.12780 -0.10640 -0.05030
## 500 meter buffer -0.10180 -0.12610 -0.01987
## 750 meter buffer -0.10070 -0.11980 -0.02426
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 4000 meter buffer 0.10910 0.210200 9 -452.803 924.4 0.00 0.174
## 3750 meter buffer 0.10790 0.194600 9 -452.890 924.6 0.17 0.160
## 3500 meter buffer 0.09872 0.174900 9 -453.403 925.6 1.20 0.096
## 4250 meter buffer 0.10780 0.214500 9 -453.501 925.8 1.39 0.087
## 3250 meter buffer 0.09166 0.172200 9 -453.536 925.9 1.47 0.084
## 250 meter buffer 0.13380 -0.007107 9 -453.650 926.1 1.69 0.075
## 4500 meter buffer 0.10720 0.214900 9 -453.722 926.3 1.84 0.069
## 4750 meter buffer 0.10540 0.211300 9 -454.084 927.0 2.56 0.048
## 3000 meter buffer 0.07873 0.163000 9 -454.247 927.3 2.89 0.041
## 2750 meter buffer 0.07848 0.157000 9 -454.350 927.5 3.09 0.037
## 5000 meter buffer 0.10130 0.211000 9 -454.427 927.7 3.25 0.034
## 2500 meter buffer 0.08318 0.131000 9 -454.921 928.7 4.24 0.021
## 1750 meter buffer 0.07713 0.098850 9 -455.199 929.2 4.79 0.016
## 2250 meter buffer 0.08659 0.119200 9 -455.334 929.5 5.06 0.014
## 2000 meter buffer 0.08297 0.105400 9 -455.336 929.5 5.07 0.014
## 1500 meter buffer 0.07840 0.082660 9 -455.732 930.3 5.86 0.009
## 1000 meter buffer 0.08467 0.066080 9 -455.971 930.8 6.34 0.007
## 1250 meter buffer 0.08330 0.063730 9 -456.055 930.9 6.50 0.007
## 500 meter buffer 0.07754 0.080650 9 -456.200 931.2 6.79 0.006
## 750 meter buffer 0.05208 0.083260 9 -457.330 933.5 9.05 0.002
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)This is quite a bit different from the global results which indicated that the 250m buffer was best fit for black bears, whereas this one it is the 4000m buffer
Top model summary
summary(bear_anthro_mods$`4000 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(black_bear, absent_black_bear) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 923.6 954.6 -452.8 905.6 223
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.003319 0.05762
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.60152 0.05278 -11.398 <2e-16 ***
## scale(harvest) 0.05542 0.06258 0.886 0.3759
## scale(seismic_lines) -0.12772 0.08033 -1.590 0.1119
## scale(seismic_lines_3D) -0.08229 0.06340 -1.298 0.1943
## scale(trails) 0.10905 0.05301 2.057 0.0397 *
## scale(wells) 0.21021 0.09903 2.123 0.0338 *
## scale(osm_industrial) -0.13807 0.06539 -2.112 0.0347 *
## scale(pipeline_transmission_lines) -0.16133 0.07616 -2.118 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Quck test on scaling
Going to run the models again WITHOUT scaling variables as recent literature suggests this may influence model results.
test_bear_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
# HFI that aren't correlated
harvest +
seismic_lines +
seismic_lines_3D +
trails +
wells +
osm_industrial +
pipeline_transmission_lines +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')model.sel(test_bear_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(hrv) cnd(osm_ind) cnd(ppl_trn_lns)
## 3750 meter buffer -0.3935 + 0.7610 -4.06200 -8.432
## 3500 meter buffer -0.4274 + 0.5648 -3.46600 -8.004
## 4250 meter buffer -0.3758 + 0.8226 -3.83500 -8.574
## 3250 meter buffer -0.4354 + 0.5984 -3.08000 -8.013
## 250 meter buffer -0.5559 + 0.0956 0.40430 -1.096
## 4500 meter buffer -0.3724 + 0.8539 -3.66900 -8.731
## 4750 meter buffer -0.3672 + 0.8794 -3.55700 -8.858
## 3000 meter buffer -0.4453 + 0.5780 -2.64100 -7.755
## 2750 meter buffer -0.4431 + 0.5598 -2.27400 -7.557
## 5000 meter buffer -0.3637 + 0.9151 -3.47500 -9.010
## 2500 meter buffer -0.4561 + 0.4954 -1.96300 -6.772
## 1750 meter buffer -0.4044 + 0.4531 -1.05300 -5.135
## 2250 meter buffer -0.4645 + 0.4344 -1.67800 -6.143
## 2000 meter buffer -0.4287 + 0.4542 -1.32400 -5.640
## 1500 meter buffer -0.4506 + 0.4292 -0.70590 -4.310
## 1000 meter buffer -0.4979 + 0.2331 -0.40810 -3.166
## 1250 meter buffer -0.4983 + 0.3435 -0.21260 -3.680
## 500 meter buffer -0.5637 + 0.2399 0.15480 -1.773
## 750 meter buffer -0.5451 + 0.2963 0.07942 -2.118
## 4000 meter buffer -0.3865 + 0.7698 -4.03600 -9.096
## cnd(ssm_lns) cnd(ssm_lns_3D) cnd(trl) cnd(wll) df logLik
## 3750 meter buffer -37.410 -7.630 80.15 12.3800 9 -452.958
## 3500 meter buffer -27.100 -7.226 71.71 9.5940 9 -453.403
## 4250 meter buffer -41.040 -7.557 82.40 12.5900 9 -453.501
## 3250 meter buffer -25.080 -7.046 65.78 9.2670 9 -453.536
## 250 meter buffer -2.274 -10.600 21.61 -0.1549 9 -453.650
## 4500 meter buffer -42.150 -7.784 83.98 12.9400 9 -453.722
## 4750 meter buffer -42.870 -8.093 84.71 13.1000 9 -454.084
## 3000 meter buffer -21.300 -7.146 55.37 8.5810 9 -454.247
## 2750 meter buffer -21.250 -7.255 53.48 8.0170 9 -454.350
## 5000 meter buffer -43.150 -8.404 83.25 13.3700 9 -454.427
## 2500 meter buffer -19.300 -7.010 54.05 6.4690 9 -454.921
## 1750 meter buffer -26.160 -7.751 41.56 4.4000 9 -455.199
## 2250 meter buffer -18.560 -6.970 54.29 5.7130 9 -455.334
## 2000 meter buffer -23.350 -7.220 48.93 4.8850 9 -455.336
## 1500 meter buffer -19.110 -8.372 37.74 3.5080 9 -455.732
## 1000 meter buffer -10.280 -9.509 33.28 2.4940 9 -455.971
## 1250 meter buffer -11.420 -9.124 36.02 2.6370 9 -456.055
## 500 meter buffer -3.477 -11.040 20.00 2.3810 9 -456.200
## 750 meter buffer -4.862 -10.420 17.61 2.8880 9 -457.330
## 4000 meter buffer -38.440 -7.300 81.98 12.9300 9
## AICc delta
## 3750 meter buffer 924.7 0.00
## 3500 meter buffer 925.6 0.89
## 4250 meter buffer 925.8 1.09
## 3250 meter buffer 925.9 1.16
## 250 meter buffer 926.1 1.38
## 4500 meter buffer 926.3 1.53
## 4750 meter buffer 927.0 2.25
## 3000 meter buffer 927.3 2.58
## 2750 meter buffer 927.5 2.78
## 5000 meter buffer 927.7 2.94
## 2500 meter buffer 928.7 3.93
## 1750 meter buffer 929.2 4.48
## 2250 meter buffer 929.5 4.75
## 2000 meter buffer 929.5 4.76
## 1500 meter buffer 930.3 5.55
## 1000 meter buffer 930.8 6.03
## 1250 meter buffer 930.9 6.19
## 500 meter buffer 931.2 6.49
## 750 meter buffer 933.5 8.74
## 4000 meter buffer
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Interesting, you get almost the same results but the 4000m buffer is the worst fit. There was a warning about model convergence issues which I suspect is the problem, this is fairly common with unscaled data. We will proceed with scaling the data as is standard practice still, but consider in discussion on ms how this may affect results.
Caribou
caribou_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(caribou, absent_caribou) ~
# HFI that aren't correlated
scale(harvest) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')# run model selection
model.sel(caribou_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(osm_ind))
## 1000 meter buffer -5.548 + -0.7434 -1.29000
## 1250 meter buffer -5.414 + -0.3827 -1.25800
## 750 meter buffer -5.861 + -1.8470 -1.09300
## 1500 meter buffer -5.316 + -0.1688 -1.08700
## 2000 meter buffer -5.303 + -0.1954 -1.32900
## 2500 meter buffer -5.357 + -0.1833 -0.98600
## 2750 meter buffer -5.357 + -0.1653 -0.65440
## 500 meter buffer -5.731 + -2.2270 -0.35600
## 1750 meter buffer -5.267 + -0.1694 -1.20600
## 2250 meter buffer -5.286 + -0.1933 -1.17600
## 3000 meter buffer -5.325 + -0.1265 -0.52740
## 3250 meter buffer -5.270 + -0.1376 -0.45680
## 3500 meter buffer -5.222 + -0.2107 -0.35290
## 250 meter buffer -205.500 + -603.1000 -0.26820
## 3750 meter buffer -5.176 + -0.3113 -0.29410
## 5000 meter buffer -5.113 + -1.1950 0.02832
## 4750 meter buffer -5.095 + -0.9945 0.01073
## 4000 meter buffer -5.124 + -0.4164 -0.21850
## 4250 meter buffer -5.081 + -0.5982 -0.13790
## 4500 meter buffer -5.069 + -0.8296 -0.06008
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 1000 meter buffer -0.3131 -0.187500 0.15460
## 1250 meter buffer -0.3097 -0.195400 0.21040
## 750 meter buffer -0.2941 -0.121300 0.22120
## 1500 meter buffer -0.4411 -0.211800 0.14590
## 2000 meter buffer -0.5624 -0.214100 0.14380
## 2500 meter buffer -0.9072 -0.206700 0.15920
## 2750 meter buffer -1.0620 -0.181200 0.10500
## 500 meter buffer -0.2392 -0.124100 0.20970
## 1750 meter buffer -0.4705 -0.196700 0.09651
## 2250 meter buffer -0.6466 -0.226500 0.17180
## 3000 meter buffer -1.1040 -0.170600 0.02474
## 3250 meter buffer -1.0460 -0.186700 -0.03486
## 3500 meter buffer -1.0620 -0.137400 -0.04571
## 250 meter buffer -0.4221 -0.179700 0.21910
## 3750 meter buffer -1.0290 -0.104100 -0.08400
## 5000 meter buffer -0.9536 0.036600 -0.03072
## 4750 meter buffer -0.9733 0.017610 -0.04811
## 4000 meter buffer -0.9482 -0.076310 -0.10210
## 4250 meter buffer -0.9276 -0.029570 -0.08694
## 4500 meter buffer -0.9202 0.003008 -0.07012
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 1000 meter buffer 0.03899 0.6921 9 -133.095 285.0 0.00 0.884
## 1250 meter buffer 0.02258 0.7364 9 -135.675 290.2 5.16 0.067
## 750 meter buffer 0.06621 0.6323 9 -136.442 291.7 6.69 0.031
## 1500 meter buffer -0.12470 0.8065 9 -138.048 294.9 9.91 0.006
## 2000 meter buffer -0.15190 0.9663 9 -139.042 296.9 11.89 0.002
## 2500 meter buffer -0.25530 1.2410 9 -139.146 297.1 12.10 0.002
## 2750 meter buffer -0.30340 1.2760 9 -139.279 297.4 12.37 0.002
## 500 meter buffer -0.03230 0.5189 9 -139.375 297.6 12.56 0.002
## 1750 meter buffer -0.13100 0.8359 9 -139.466 297.7 12.74 0.002
## 2250 meter buffer -0.17320 1.0390 9 -139.996 298.8 13.80 0.001
## 3000 meter buffer -0.39110 1.2520 9 -140.152 299.1 14.12 0.001
## 3250 meter buffer -0.44820 1.1740 9 -141.535 301.9 16.88 0.000
## 3500 meter buffer -0.48140 1.1340 9 -142.032 302.9 17.87 0.000
## 250 meter buffer -0.26720 0.3502 9 -142.102 303.0 18.01 0.000
## 3750 meter buffer -0.45830 1.0850 9 -142.635 304.1 19.08 0.000
## 5000 meter buffer -0.54130 0.8192 9 -143.734 306.3 21.28 0.000
## 4750 meter buffer -0.53940 0.8414 9 -143.859 306.5 21.53 0.000
## 4000 meter buffer -0.40610 0.9949 9 -143.916 306.6 21.64 0.000
## 4250 meter buffer -0.42730 0.9081 9 -144.526 307.9 22.86 0.000
## 4500 meter buffer -0.45820 0.8304 9 -144.675 308.2 23.16 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)Let’s look at the summary for the top model
summary(caribou_anthro_mods$`1000 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(caribou, absent_caribou) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 284.2 315.2 -133.1 266.2 223
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 4.648 2.156
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.54820 1.05969 -5.236 1.64e-07 ***
## scale(harvest) -0.74343 0.66903 -1.111 0.2665
## scale(seismic_lines) 0.15460 0.16295 0.949 0.3428
## scale(seismic_lines_3D) -0.18749 0.19623 -0.955 0.3393
## scale(trails) 0.03899 0.20948 0.186 0.8524
## scale(wells) 0.69213 0.15941 4.342 1.41e-05 ***
## scale(osm_industrial) -1.28957 0.45284 -2.848 0.0044 **
## scale(pipeline_transmission_lines) -0.31311 0.28058 -1.116 0.2645
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Coyote
coyote_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(coyote, absent_coyote) ~
# HFI that aren't correlated
scale(harvest) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(coyote_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(osm_ind))
## 4750 meter buffer -1.393 + -0.168500 0.4460
## 5000 meter buffer -1.394 + -0.178600 0.4404
## 4500 meter buffer -1.391 + -0.149400 0.4498
## 3000 meter buffer -1.382 + -0.073710 0.4564
## 2750 meter buffer -1.380 + -0.067890 0.4562
## 4250 meter buffer -1.388 + -0.132400 0.4489
## 3250 meter buffer -1.382 + -0.082540 0.4497
## 2250 meter buffer -1.379 + -0.049800 0.4119
## 2500 meter buffer -1.379 + -0.057110 0.4266
## 4000 meter buffer -1.386 + -0.124100 0.4384
## 2000 meter buffer -1.378 + -0.034050 0.4143
## 1750 meter buffer -1.379 + -0.017100 0.4133
## 3500 meter buffer -1.383 + -0.099850 0.4312
## 3750 meter buffer -1.384 + -0.110800 0.4269
## 1500 meter buffer -1.381 + -0.004780 0.4114
## 1250 meter buffer -1.383 + -0.009980 0.3993
## 750 meter buffer -1.388 + 0.019510 0.3722
## 1000 meter buffer -1.385 + -0.013770 0.3735
## 500 meter buffer -1.388 + 0.008353 0.3456
## 250 meter buffer -1.377 + -0.091720 0.2663
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 4750 meter buffer 0.0492900 -0.05031 0.1597
## 5000 meter buffer 0.0554700 -0.04634 0.1629
## 4500 meter buffer 0.0288900 -0.04988 0.1558
## 3000 meter buffer -0.0458300 -0.05317 0.1214
## 2750 meter buffer -0.0361000 -0.06249 0.1268
## 4250 meter buffer 0.0123000 -0.05042 0.1396
## 3250 meter buffer -0.0367900 -0.04980 0.1168
## 2250 meter buffer 0.0258900 -0.10060 0.1738
## 2500 meter buffer 0.0004534 -0.07837 0.1647
## 4000 meter buffer -0.0152200 -0.04670 0.1304
## 2000 meter buffer 0.0341400 -0.12410 0.1496
## 1750 meter buffer 0.0524900 -0.14990 0.1417
## 3500 meter buffer -0.0254300 -0.04930 0.1189
## 3750 meter buffer -0.0200100 -0.04823 0.1191
## 1500 meter buffer 0.0590100 -0.17360 0.1518
## 1250 meter buffer 0.0761600 -0.18780 0.1577
## 750 meter buffer 0.1324000 -0.15780 0.1663
## 1000 meter buffer 0.1118000 -0.16520 0.1550
## 500 meter buffer 0.0713100 -0.12700 0.1541
## 250 meter buffer 0.1284000 -0.05137 0.1034
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 4750 meter buffer 0.021530 0.15120 9 -478.831 976.5 0.00 0.268
## 5000 meter buffer 0.030250 0.15510 9 -478.960 976.7 0.26 0.235
## 4500 meter buffer 0.011280 0.16170 9 -479.384 977.6 1.11 0.154
## 3000 meter buffer 0.020230 0.22810 9 -479.674 978.2 1.69 0.115
## 2750 meter buffer 0.018040 0.22270 9 -480.403 979.6 3.15 0.056
## 4250 meter buffer 0.004725 0.16960 9 -480.460 979.7 3.26 0.052
## 3250 meter buffer 0.013730 0.21790 9 -480.691 980.2 3.72 0.042
## 2250 meter buffer 0.058870 0.24090 9 -481.638 982.1 5.61 0.016
## 2500 meter buffer 0.040200 0.23510 9 -481.652 982.1 5.64 0.016
## 4000 meter buffer 0.009785 0.19210 9 -481.857 982.5 6.05 0.013
## 2000 meter buffer 0.050720 0.23600 9 -482.070 983.0 6.48 0.010
## 1750 meter buffer 0.049730 0.23850 9 -482.384 983.6 7.11 0.008
## 3500 meter buffer 0.021720 0.21730 9 -482.448 983.7 7.23 0.007
## 3750 meter buffer 0.020430 0.21150 9 -482.925 984.7 8.19 0.004
## 1500 meter buffer 0.038900 0.25470 9 -483.202 985.2 8.74 0.003
## 1250 meter buffer 0.018110 0.23530 9 -487.827 994.5 17.99 0.000
## 750 meter buffer -0.073220 0.28850 9 -493.605 1006.0 29.55 0.000
## 1000 meter buffer 0.002133 0.24340 9 -493.978 1006.8 30.29 0.000
## 500 meter buffer -0.113300 0.24890 9 -499.558 1017.9 41.45 0.000
## 250 meter buffer 0.008134 0.03331 9 -516.501 1051.8 75.34 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(coyote_anthro_mods$`4750 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(coyote, absent_coyote) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 975.7 1006.7 -478.8 957.7 223
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.3289 0.5735
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.39250 0.24038 -5.793 6.92e-09 ***
## scale(harvest) -0.16846 0.09498 -1.774 0.0761 .
## scale(seismic_lines) 0.15970 0.07739 2.064 0.0391 *
## scale(seismic_lines_3D) -0.05031 0.07036 -0.715 0.4746
## scale(trails) 0.02153 0.06545 0.329 0.7421
## scale(wells) 0.15121 0.08582 1.762 0.0781 .
## scale(osm_industrial) 0.44601 0.06039 7.386 1.51e-13 ***
## scale(pipeline_transmission_lines) 0.04929 0.08745 0.564 0.5730
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fisher
fisher_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(fisher, absent_fisher) ~
# HFI that aren't correlated
scale(harvest) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(fisher_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(osm_ind))
## 250 meter buffer -2.970 + -0.37280 -0.261500
## 500 meter buffer -2.969 + -0.32480 -0.235700
## 2500 meter buffer -2.951 + -0.06196 0.058010
## 750 meter buffer -2.954 + -0.25240 -0.099590
## 2750 meter buffer -2.948 + -0.05096 0.075520
## 3000 meter buffer -2.943 + -0.04421 0.080270
## 3250 meter buffer -2.940 + -0.06051 0.044190
## 1000 meter buffer -2.956 + -0.24610 -0.034300
## 1250 meter buffer -2.960 + -0.26160 -0.018860
## 4500 meter buffer -2.953 + -0.14600 -0.238000
## 4750 meter buffer -2.958 + -0.17650 -0.291000
## 5000 meter buffer -2.961 + -0.20230 -0.321400
## 4250 meter buffer -2.948 + -0.12070 -0.176100
## 4000 meter buffer -2.942 + -0.10750 -0.085540
## 2250 meter buffer -2.946 + -0.09139 0.066760
## 3500 meter buffer -2.938 + -0.08978 0.002052
## 3750 meter buffer -2.940 + -0.10760 -0.040320
## 2000 meter buffer -2.944 + -0.12840 0.059920
## 1500 meter buffer -2.951 + -0.21610 0.028080
## 1750 meter buffer -2.947 + -0.17880 0.063700
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 250 meter buffer -0.1094 -0.1564 -0.05936
## 500 meter buffer -0.2236 -0.2159 -0.05955
## 2500 meter buffer -0.3526 -0.2064 -0.12610
## 750 meter buffer -0.2420 -0.2620 -0.02603
## 2750 meter buffer -0.3390 -0.1925 -0.15160
## 3000 meter buffer -0.3117 -0.1849 -0.14380
## 3250 meter buffer -0.2765 -0.1895 -0.15020
## 1000 meter buffer -0.2149 -0.3244 -0.03179
## 1250 meter buffer -0.1835 -0.3590 -0.07089
## 4500 meter buffer -0.1914 -0.2183 -0.12830
## 4750 meter buffer -0.1599 -0.2316 -0.11390
## 5000 meter buffer -0.1183 -0.2550 -0.11200
## 4250 meter buffer -0.2046 -0.2154 -0.13470
## 4000 meter buffer -0.2241 -0.2092 -0.12610
## 2250 meter buffer -0.3131 -0.2203 -0.09885
## 3500 meter buffer -0.2190 -0.1963 -0.15260
## 3750 meter buffer -0.2047 -0.2082 -0.14090
## 2000 meter buffer -0.2581 -0.2525 -0.10550
## 1500 meter buffer -0.1904 -0.3317 -0.07759
## 1750 meter buffer -0.2026 -0.2945 -0.09132
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 250 meter buffer -0.004209 -0.1091 9 -282.462 583.7 0.00 0.395
## 500 meter buffer 0.004649 0.1249 9 -282.654 584.1 0.38 0.326
## 2500 meter buffer -0.027110 0.4010 9 -284.705 588.2 4.49 0.042
## 750 meter buffer -0.009610 0.1642 9 -284.903 588.6 4.88 0.034
## 2750 meter buffer -0.048780 0.3795 9 -285.001 588.8 5.08 0.031
## 3000 meter buffer -0.069300 0.3749 9 -285.410 589.6 5.90 0.021
## 3250 meter buffer -0.078650 0.3791 9 -285.413 589.6 5.90 0.021
## 1000 meter buffer 0.078410 0.1735 9 -285.595 590.0 6.27 0.017
## 1250 meter buffer 0.078080 0.2060 9 -285.769 590.3 6.61 0.014
## 4500 meter buffer 0.026860 0.4652 9 -285.857 590.5 6.79 0.013
## 4750 meter buffer 0.063950 0.4717 9 -285.860 590.5 6.80 0.013
## 5000 meter buffer 0.084000 0.4533 9 -285.905 590.6 6.89 0.013
## 4250 meter buffer -0.000520 0.4439 9 -285.963 590.7 7.00 0.012
## 4000 meter buffer -0.008931 0.4188 9 -286.136 591.1 7.35 0.010
## 2250 meter buffer -0.033650 0.3563 9 -286.186 591.2 7.45 0.010
## 3500 meter buffer -0.033270 0.3558 9 -286.281 591.4 7.64 0.009
## 3750 meter buffer -0.004305 0.3783 9 -286.283 591.4 7.64 0.009
## 2000 meter buffer -0.023340 0.3050 9 -286.935 592.7 8.95 0.005
## 1500 meter buffer 0.061120 0.2195 9 -287.165 593.1 9.41 0.004
## 1750 meter buffer 0.042090 0.2366 9 -287.421 593.7 9.92 0.003
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(fisher_anthro_mods$`250 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(fisher, absent_fisher) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 582.9 613.9 -282.5 564.9 223
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.7269 0.8526
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.970098 0.365121 -8.135 4.13e-16 ***
## scale(harvest) -0.372810 0.115254 -3.235 0.00122 **
## scale(seismic_lines) -0.059365 0.081191 -0.731 0.46467
## scale(seismic_lines_3D) -0.156392 0.127570 -1.226 0.22023
## scale(trails) -0.004209 0.091425 -0.046 0.96328
## scale(wells) -0.109099 0.081448 -1.339 0.18041
## scale(osm_industrial) -0.261515 0.119236 -2.193 0.02829 *
## scale(pipeline_transmission_lines) -0.109381 0.130572 -0.838 0.40220
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Grey wolf
wolf_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(grey_wolf, absent_grey_wolf) ~
# HFI that aren't correlated
scale(harvest) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(wolf_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(osm_ind))
## 2000 meter buffer -3.221 + 0.10340 -0.32840
## 2250 meter buffer -3.226 + 0.12440 -0.29220
## 1750 meter buffer -3.211 + 0.08444 -0.29030
## 2500 meter buffer -3.223 + 0.14350 -0.23780
## 1500 meter buffer -3.194 + 0.07258 -0.27170
## 2750 meter buffer -3.214 + 0.16500 -0.19800
## 3000 meter buffer -3.202 + 0.16030 -0.19160
## 1250 meter buffer -3.176 + 0.07269 -0.24850
## 3250 meter buffer -3.194 + 0.15380 -0.18540
## 3500 meter buffer -3.182 + 0.14850 -0.18730
## 3750 meter buffer -3.175 + 0.13860 -0.17400
## 1000 meter buffer -3.168 + 0.06317 -0.18990
## 4000 meter buffer -3.168 + 0.13200 -0.15770
## 4250 meter buffer -3.162 + 0.10390 -0.13340
## 4500 meter buffer -3.156 + 0.08534 -0.11930
## 4750 meter buffer -3.150 + 0.06932 -0.11030
## 500 meter buffer -3.155 + -0.00590 -0.04929
## 250 meter buffer -3.197 + 0.03152 0.05259
## 5000 meter buffer -3.145 + 0.05273 -0.10220
## 750 meter buffer -3.146 + 0.07042 -0.11810
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 2000 meter buffer 0.115100 -0.8980 -0.02981
## 2250 meter buffer 0.116300 -0.9238 -0.02200
## 1750 meter buffer 0.086980 -0.8791 -0.10100
## 2500 meter buffer 0.118500 -0.9636 -0.01864
## 1500 meter buffer 0.051250 -0.8376 -0.06989
## 2750 meter buffer 0.098000 -0.9549 -0.02822
## 3000 meter buffer 0.092660 -0.9213 -0.03185
## 1250 meter buffer 0.029690 -0.7955 -0.01897
## 3250 meter buffer 0.090920 -0.8858 -0.04028
## 3500 meter buffer 0.075830 -0.8552 -0.04462
## 3750 meter buffer 0.046680 -0.8153 -0.02611
## 1000 meter buffer 0.032470 -0.8060 0.03936
## 4000 meter buffer -0.010080 -0.7769 -0.02803
## 4250 meter buffer -0.047900 -0.7534 -0.02356
## 4500 meter buffer -0.080920 -0.7363 -0.03010
## 4750 meter buffer -0.120000 -0.7114 -0.02722
## 500 meter buffer -0.026250 -0.8061 0.06533
## 250 meter buffer -0.005518 -1.0380 0.01478
## 5000 meter buffer -0.153000 -0.6968 -0.02332
## 750 meter buffer 0.051550 -0.8092 0.11840
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 2000 meter buffer 0.22940 -0.19410 9 -245.933 510.7 0.00 0.245
## 2250 meter buffer 0.21720 -0.22140 9 -245.978 510.8 0.09 0.234
## 1750 meter buffer 0.21200 -0.18940 9 -246.342 511.5 0.82 0.163
## 2500 meter buffer 0.20780 -0.20800 9 -246.597 512.0 1.33 0.126
## 1500 meter buffer 0.19670 -0.17190 9 -247.220 513.3 2.57 0.068
## 2750 meter buffer 0.17720 -0.20480 9 -247.504 513.8 3.14 0.051
## 3000 meter buffer 0.16120 -0.21970 9 -248.095 515.0 4.32 0.028
## 1250 meter buffer 0.15410 -0.19940 9 -248.495 515.8 5.12 0.019
## 3250 meter buffer 0.15590 -0.25630 9 -248.532 515.9 5.20 0.018
## 3500 meter buffer 0.13720 -0.25200 9 -249.144 517.1 6.42 0.010
## 3750 meter buffer 0.11180 -0.27550 9 -249.467 517.7 7.07 0.007
## 1000 meter buffer 0.10150 -0.25290 9 -249.553 517.9 7.24 0.007
## 4000 meter buffer 0.08157 -0.25910 9 -249.716 518.2 7.57 0.006
## 4250 meter buffer 0.07508 -0.26330 9 -249.998 518.8 8.13 0.004
## 4500 meter buffer 0.06654 -0.24130 9 -250.309 519.4 8.75 0.003
## 4750 meter buffer 0.05308 -0.22110 9 -250.448 519.7 9.03 0.003
## 500 meter buffer 0.10580 -0.32170 9 -250.493 519.8 9.12 0.003
## 250 meter buffer 0.22240 -0.09401 9 -250.544 519.9 9.22 0.002
## 5000 meter buffer 0.04544 -0.19760 9 -250.565 519.9 9.26 0.002
## 750 meter buffer 0.03785 -0.14140 9 -251.677 522.2 11.49 0.001
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(wolf_anthro_mods$`2000 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(grey_wolf, absent_grey_wolf) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 509.9 540.9 -245.9 491.9 223
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.1494 0.3865
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.22060 0.20918 -15.396 < 2e-16 ***
## scale(harvest) 0.10337 0.11432 0.904 0.36588
## scale(seismic_lines) -0.02981 0.11606 -0.257 0.79729
## scale(seismic_lines_3D) -0.89796 0.29493 -3.045 0.00233 **
## scale(trails) 0.22937 0.10300 2.227 0.02595 *
## scale(wells) -0.19412 0.19809 -0.980 0.32711
## scale(osm_industrial) -0.32843 0.18340 -1.791 0.07333 .
## scale(pipeline_transmission_lines) 0.11514 0.13458 0.856 0.39225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Lynx
lynx_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(lynx, absent_lynx) ~
# HFI that aren't correlated
scale(harvest) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(lynx_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(osm_ind))
## 250 meter buffer -1.895 + -0.29470 0.09425
## 500 meter buffer -1.873 + -0.21570 0.10290
## 4250 meter buffer -1.878 + 0.17710 0.05193
## 5000 meter buffer -1.875 + 0.13110 0.08080
## 4750 meter buffer -1.875 + 0.14570 0.07245
## 3750 meter buffer -1.878 + 0.16610 0.04198
## 4500 meter buffer -1.877 + 0.17010 0.06234
## 3500 meter buffer -1.876 + 0.15980 0.04248
## 4000 meter buffer -1.878 + 0.17160 0.04480
## 1000 meter buffer -1.866 + -0.11440 0.07094
## 3250 meter buffer -1.875 + 0.15400 0.04847
## 750 meter buffer -1.867 + -0.15050 0.09380
## 3000 meter buffer -1.874 + 0.14250 0.04267
## 1250 meter buffer -1.866 + -0.10470 0.05542
## 2750 meter buffer -1.871 + 0.12220 0.04221
## 1500 meter buffer -1.864 + -0.07129 0.04894
## 2000 meter buffer -1.865 + 0.03122 0.04712
## 2500 meter buffer -1.868 + 0.09171 0.04379
## 1750 meter buffer -1.863 + -0.01427 0.04891
## 2250 meter buffer -1.866 + 0.06938 0.04817
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 250 meter buffer 0.0066170 0.07458 -0.165700
## 500 meter buffer -0.0196700 0.06938 -0.081080
## 4250 meter buffer -0.0350500 0.07413 -0.190600
## 5000 meter buffer -0.0580300 0.10100 -0.175100
## 4750 meter buffer -0.0504500 0.09102 -0.172100
## 3750 meter buffer -0.0222900 0.06483 -0.195500
## 4500 meter buffer -0.0506400 0.08297 -0.178900
## 3500 meter buffer -0.0123700 0.06459 -0.181900
## 4000 meter buffer -0.0370000 0.06961 -0.190900
## 1000 meter buffer 0.0197900 0.07077 -0.006769
## 3250 meter buffer -0.0129000 0.07140 -0.166700
## 750 meter buffer 0.0026100 0.06941 -0.037170
## 3000 meter buffer -0.0126300 0.07540 -0.145000
## 1250 meter buffer 0.0222000 0.07157 0.039080
## 2750 meter buffer -0.0067840 0.07739 -0.123700
## 1500 meter buffer 0.0266300 0.07479 0.036140
## 2000 meter buffer 0.0283600 0.08070 -0.001350
## 2500 meter buffer -0.0073830 0.07901 -0.097930
## 1750 meter buffer 0.0282200 0.07890 0.016560
## 2250 meter buffer 0.0007309 0.08147 -0.063070
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 250 meter buffer 0.06159 -0.2237 9 -437.719 894.2 0.00 0.996
## 500 meter buffer -0.01431 -0.2005 9 -444.187 907.2 12.94 0.002
## 4250 meter buffer 0.09287 -0.1310 9 -445.767 910.3 16.10 0.000
## 5000 meter buffer 0.08851 -0.1159 9 -445.873 910.6 16.31 0.000
## 4750 meter buffer 0.09413 -0.1196 9 -445.948 910.7 16.46 0.000
## 3750 meter buffer 0.09472 -0.1410 9 -445.969 910.7 16.50 0.000
## 4500 meter buffer 0.09220 -0.1131 9 -446.011 910.8 16.58 0.000
## 3500 meter buffer 0.09696 -0.1499 9 -446.187 911.2 16.94 0.000
## 4000 meter buffer 0.09032 -0.1231 9 -446.224 911.3 17.01 0.000
## 1000 meter buffer 0.09300 -0.2122 9 -446.379 911.6 17.32 0.000
## 3250 meter buffer 0.08939 -0.1602 9 -446.521 911.9 17.60 0.000
## 750 meter buffer 0.03048 -0.1860 9 -446.922 912.7 18.41 0.000
## 3000 meter buffer 0.08870 -0.1684 9 -447.049 912.9 18.66 0.000
## 1250 meter buffer 0.10150 -0.2185 9 -447.156 913.1 18.87 0.000
## 2750 meter buffer 0.07960 -0.1789 9 -447.989 914.8 20.54 0.000
## 1500 meter buffer 0.10920 -0.2021 9 -448.355 915.5 21.27 0.000
## 2000 meter buffer 0.10340 -0.2143 9 -448.722 916.3 22.01 0.000
## 2500 meter buffer 0.08089 -0.1775 9 -448.757 916.3 22.08 0.000
## 1750 meter buffer 0.10430 -0.2096 9 -448.838 916.5 22.24 0.000
## 2250 meter buffer 0.08931 -0.1882 9 -448.923 916.7 22.41 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(lynx_anthro_mods$`250 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(lynx, absent_lynx) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 893.4 924.5 -437.7 875.4 223
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.1606 0.4007
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.894836 0.174410 -10.864 < 2e-16 ***
## scale(harvest) -0.294671 0.090226 -3.266 0.00109 **
## scale(seismic_lines) -0.165728 0.071320 -2.324 0.02014 *
## scale(seismic_lines_3D) 0.074577 0.053207 1.402 0.16102
## scale(trails) 0.061593 0.052606 1.171 0.24166
## scale(wells) -0.223700 0.073888 -3.028 0.00247 **
## scale(osm_industrial) 0.094250 0.051713 1.823 0.06837 .
## scale(pipeline_transmission_lines) 0.006617 0.059660 0.111 0.91169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Moose
moose_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(moose, absent_moose) ~
# HFI that aren't correlated
scale(harvest) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(moose_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(osm_ind))
## 750 meter buffer -1.835 + 0.0183700 -0.3434
## 1000 meter buffer -1.831 + -0.0007321 -0.3139
## 500 meter buffer -1.827 + 0.0069890 -0.2877
## 1250 meter buffer -1.822 + 0.0138900 -0.2539
## 1500 meter buffer -1.822 + 0.0004848 -0.2435
## 2500 meter buffer -1.819 + -0.0374500 -0.2053
## 3000 meter buffer -1.818 + -0.0412000 -0.2303
## 1750 meter buffer -1.819 + -0.0147400 -0.1903
## 2750 meter buffer -1.818 + -0.0409500 -0.2233
## 3250 meter buffer -1.817 + -0.0479700 -0.2407
## 3500 meter buffer -1.816 + -0.0602400 -0.2486
## 3750 meter buffer -1.815 + -0.0718000 -0.2567
## 2000 meter buffer -1.818 + -0.0282200 -0.1659
## 2250 meter buffer -1.818 + -0.0252900 -0.1770
## 4000 meter buffer -1.814 + -0.0782900 -0.2503
## 250 meter buffer -1.814 + -0.0025580 -0.1621
## 4250 meter buffer -1.812 + -0.0881700 -0.2331
## 4500 meter buffer -1.811 + -0.0990100 -0.2128
## 5000 meter buffer -1.810 + -0.1124000 -0.1833
## 4750 meter buffer -1.810 + -0.1050000 -0.1909
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 750 meter buffer 0.0966400 0.087410 0.015410
## 1000 meter buffer 0.0793700 0.068490 -0.050990
## 500 meter buffer 0.0700700 0.100300 0.069440
## 1250 meter buffer 0.0877400 0.054910 -0.038170
## 1500 meter buffer 0.0886300 0.045150 -0.065370
## 2500 meter buffer 0.1786000 -0.024930 -0.012030
## 3000 meter buffer 0.1953000 -0.032250 -0.030410
## 1750 meter buffer 0.1020000 0.027390 -0.075250
## 2750 meter buffer 0.1921000 -0.032160 -0.014060
## 3250 meter buffer 0.1915000 -0.033870 -0.030790
## 3500 meter buffer 0.1885000 -0.035860 -0.029600
## 3750 meter buffer 0.1807000 -0.037740 -0.026340
## 2000 meter buffer 0.1176000 0.008666 -0.042810
## 2250 meter buffer 0.1450000 -0.008776 -0.022820
## 4000 meter buffer 0.1638000 -0.037310 -0.026490
## 250 meter buffer -0.0008117 0.111900 0.090400
## 4250 meter buffer 0.1441000 -0.032050 -0.013890
## 4500 meter buffer 0.1357000 -0.029150 0.002179
## 5000 meter buffer 0.1556000 -0.033120 0.019220
## 4750 meter buffer 0.1413000 -0.028300 0.014930
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 750 meter buffer 0.05467 -0.15570 9 -433.551 885.9 0.00 0.526
## 1000 meter buffer 0.07086 -0.20550 9 -434.001 886.8 0.90 0.335
## 500 meter buffer 0.06798 -0.09246 9 -435.787 890.4 4.47 0.056
## 1250 meter buffer 0.08575 -0.18130 9 -437.308 893.4 7.51 0.012
## 1500 meter buffer 0.09491 -0.20160 9 -437.362 893.5 7.62 0.012
## 2500 meter buffer 0.09335 -0.25700 9 -437.712 894.2 8.32 0.008
## 3000 meter buffer 0.09252 -0.22390 9 -437.745 894.3 8.39 0.008
## 1750 meter buffer 0.07996 -0.25640 9 -437.870 894.6 8.64 0.007
## 2750 meter buffer 0.09193 -0.22680 9 -437.898 894.6 8.69 0.007
## 3250 meter buffer 0.09264 -0.19940 9 -438.052 894.9 9.00 0.006
## 3500 meter buffer 0.09262 -0.18520 9 -438.227 895.3 9.35 0.005
## 3750 meter buffer 0.09553 -0.15720 9 -438.516 895.8 9.93 0.004
## 2000 meter buffer 0.07710 -0.27730 9 -438.578 896.0 10.05 0.003
## 2250 meter buffer 0.08022 -0.26710 9 -438.593 896.0 10.08 0.003
## 4000 meter buffer 0.09752 -0.14630 9 -438.875 896.6 10.65 0.003
## 250 meter buffer 0.09090 -0.03078 9 -439.130 897.1 11.16 0.002
## 4250 meter buffer 0.10080 -0.14190 9 -439.641 898.1 12.18 0.001
## 4500 meter buffer 0.10270 -0.14890 9 -440.156 899.1 13.21 0.001
## 5000 meter buffer 0.10460 -0.18440 9 -440.256 899.3 13.41 0.001
## 4750 meter buffer 0.10210 -0.17380 9 -440.331 899.5 13.56 0.001
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(moose_anthro_mods$`750 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(moose, absent_moose) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 885.1 916.1 -433.6 867.1 223
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.2396 0.4895
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.83485 0.20929 -8.767 < 2e-16 ***
## scale(harvest) 0.01837 0.06765 0.272 0.785981
## scale(seismic_lines) 0.01541 0.06597 0.234 0.815250
## scale(seismic_lines_3D) 0.08741 0.06029 1.450 0.147104
## scale(trails) 0.05467 0.05245 1.042 0.297238
## scale(wells) -0.15573 0.08944 -1.741 0.081644 .
## scale(osm_industrial) -0.34339 0.08957 -3.834 0.000126 ***
## scale(pipeline_transmission_lines) 0.09664 0.06979 1.385 0.166093
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Red fox
fox_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(red_fox, absent_red_fox) ~
# HFI that aren't correlated
scale(harvest) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(fox_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(osm_ind))
## 1750 meter buffer -4.593 + -0.26290 0.44700
## 2000 meter buffer -4.533 + -0.21720 0.39520
## 1500 meter buffer -4.486 + -0.37670 0.36810
## 2250 meter buffer -4.490 + -0.24490 0.36720
## 2500 meter buffer -4.509 + -0.19350 0.39310
## 4250 meter buffer -4.480 + -0.13870 0.28240
## 4500 meter buffer -4.514 + -0.14030 0.42600
## 1250 meter buffer -4.370 + -0.39910 0.31620
## 4000 meter buffer -4.453 + -0.15280 0.11610
## 3750 meter buffer -4.457 + -0.16130 0.07543
## 4750 meter buffer -4.511 + -0.12340 0.44770
## 2750 meter buffer -4.467 + -0.14090 0.35290
## 5000 meter buffer -4.469 + -0.08203 0.42310
## 3000 meter buffer -4.452 + -0.12980 0.30140
## 500 meter buffer -4.374 + -0.18810 0.38670
## 3250 meter buffer -4.408 + -0.14980 0.24890
## 3500 meter buffer -4.427 + -0.17300 0.21850
## 1000 meter buffer -4.310 + -0.38190 0.19830
## 750 meter buffer -4.324 + -0.37680 0.30650
## 250 meter buffer -4.307 + -0.15800 0.24810
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 1750 meter buffer 0.324000 -0.22090 -0.4869
## 2000 meter buffer 0.345600 -0.29500 -0.4768
## 1500 meter buffer 0.217900 -0.12130 -0.4365
## 2250 meter buffer 0.403800 -0.33790 -0.3919
## 2500 meter buffer 0.466100 -0.36660 -0.4114
## 4250 meter buffer 0.714900 -0.69650 -0.2710
## 4500 meter buffer 0.753200 -0.62740 -0.3167
## 1250 meter buffer 0.059680 -0.13260 -0.3583
## 4000 meter buffer 0.719400 -0.67060 -0.2413
## 3750 meter buffer 0.758500 -0.63570 -0.2504
## 4750 meter buffer 0.766800 -0.65750 -0.3209
## 2750 meter buffer 0.502500 -0.41890 -0.3420
## 5000 meter buffer 0.655300 -0.71140 -0.2513
## 3000 meter buffer 0.600900 -0.42730 -0.2605
## 500 meter buffer 0.054000 -0.01504 -0.2116
## 3250 meter buffer 0.562200 -0.47190 -0.1988
## 3500 meter buffer 0.702400 -0.50120 -0.2700
## 1000 meter buffer 0.008758 -0.09754 -0.3227
## 750 meter buffer 0.130300 -0.07750 -0.1869
## 250 meter buffer -0.033860 0.12030 -0.1493
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 1750 meter buffer -0.8835 -0.08828 9 -173.699 366.2 0.00 0.733
## 2000 meter buffer -0.7528 0.01459 9 -175.241 369.3 3.08 0.157
## 1500 meter buffer -0.7755 -0.03183 9 -176.233 371.3 5.07 0.058
## 2250 meter buffer -0.6416 0.04205 9 -177.849 374.5 8.30 0.012
## 2500 meter buffer -0.6323 0.01523 9 -178.330 375.5 9.26 0.007
## 4250 meter buffer -0.3811 0.21500 9 -178.624 376.1 9.85 0.005
## 4500 meter buffer -0.4312 0.09273 9 -178.786 376.4 10.17 0.005
## 1250 meter buffer -0.6934 0.19530 9 -178.788 376.4 10.18 0.005
## 4000 meter buffer -0.2684 0.27960 9 -178.884 376.6 10.37 0.004
## 3750 meter buffer -0.1880 0.25870 9 -178.933 376.7 10.47 0.004
## 4750 meter buffer -0.3783 0.08301 9 -179.146 377.1 10.89 0.003
## 2750 meter buffer -0.4960 0.11540 9 -179.253 377.3 11.11 0.003
## 5000 meter buffer -0.3243 0.21080 9 -179.989 378.8 12.58 0.001
## 3000 meter buffer -0.3654 0.16090 9 -180.057 378.9 12.72 0.001
## 500 meter buffer -0.6258 0.47560 9 -180.174 379.2 12.95 0.001
## 3250 meter buffer -0.2426 0.25470 9 -180.556 379.9 13.71 0.001
## 3500 meter buffer -0.1972 0.14450 9 -180.576 380.0 13.75 0.001
## 1000 meter buffer -0.6295 0.25200 9 -182.729 384.3 18.06 0.000
## 750 meter buffer -0.6376 0.29360 9 -183.234 385.3 19.07 0.000
## 250 meter buffer -0.8490 0.13910 9 -186.735 392.3 26.07 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(fox_anthro_mods$`1750 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(red_fox, absent_red_fox) ~ scale(harvest) + scale(seismic_lines) +
## scale(seismic_lines_3D) + scale(trails) + scale(wells) +
## scale(osm_industrial) + scale(pipeline_transmission_lines) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 365.4 396.4 -173.7 347.4 223
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 2.963 1.721
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.59261 0.76036 -6.040 1.54e-09 ***
## scale(harvest) -0.26287 0.24399 -1.077 0.281304
## scale(seismic_lines) -0.48687 0.15878 -3.066 0.002167 **
## scale(seismic_lines_3D) -0.22086 0.29037 -0.761 0.446900
## scale(trails) -0.88347 0.25707 -3.437 0.000589 ***
## scale(wells) -0.08828 0.24232 -0.364 0.715610
## scale(osm_industrial) 0.44701 0.17367 2.574 0.010054 *
## scale(pipeline_transmission_lines) 0.32404 0.23535 1.377 0.168557
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1White-tailed deer
deer_anthro_mods <- osm_anthro_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~
# HFI that aren't correlated
scale(harvest) +
scale(seismic_lines) +
scale(seismic_lines_3D) +
scale(trails) +
scale(wells) +
scale(osm_industrial) +
scale(pipeline_transmission_lines) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(deer_anthro_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(osm_ind))
## 1500 meter buffer -0.1950 + -0.08754 0.3267
## 1250 meter buffer -0.1941 + -0.08803 0.3532
## 1750 meter buffer -0.1896 + -0.10390 0.2888
## 4000 meter buffer -0.1971 + -0.23130 0.2007
## 3750 meter buffer -0.1931 + -0.20820 0.1936
## 2000 meter buffer -0.1887 + -0.11210 0.2603
## 2250 meter buffer -0.1851 + -0.12970 0.2433
## 4250 meter buffer -0.2004 + -0.24410 0.2087
## 3500 meter buffer -0.1921 + -0.18820 0.1937
## 3250 meter buffer -0.1879 + -0.16530 0.1997
## 4500 meter buffer -0.2040 + -0.25290 0.2178
## 2500 meter buffer -0.1853 + -0.13080 0.2278
## 2750 meter buffer -0.1840 + -0.13380 0.2179
## 1000 meter buffer -0.1872 + -0.11660 0.3465
## 3000 meter buffer -0.1846 + -0.15090 0.2114
## 4750 meter buffer -0.2082 + -0.27160 0.2229
## 5000 meter buffer -0.2116 + -0.29100 0.2223
## 750 meter buffer -0.1785 + -0.08545 0.3466
## 500 meter buffer -0.1836 + -0.01532 0.3324
## 250 meter buffer -0.1848 + 0.02891 0.2452
## cnd(scl(ppl_trn_lns)) cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns))
## 1500 meter buffer 0.068920 -0.4287 -0.07579
## 1250 meter buffer 0.071360 -0.4240 -0.06643
## 1750 meter buffer 0.073990 -0.4028 -0.09030
## 4000 meter buffer 0.016820 -0.3929 -0.03263
## 3750 meter buffer 0.046180 -0.4028 -0.02544
## 2000 meter buffer 0.058620 -0.3906 -0.08964
## 2250 meter buffer 0.048600 -0.3651 -0.05133
## 4250 meter buffer -0.003816 -0.3771 -0.05001
## 3500 meter buffer 0.047900 -0.3998 -0.03284
## 3250 meter buffer 0.055870 -0.3867 -0.02725
## 4500 meter buffer -0.019730 -0.3706 -0.04690
## 2500 meter buffer 0.048840 -0.3605 -0.05163
## 2750 meter buffer 0.045340 -0.3666 -0.03939
## 1000 meter buffer 0.091750 -0.3536 -0.07067
## 3000 meter buffer 0.050280 -0.3692 -0.02601
## 4750 meter buffer -0.046200 -0.3562 -0.03454
## 5000 meter buffer -0.087870 -0.3352 -0.02876
## 750 meter buffer 0.127200 -0.2965 -0.02425
## 500 meter buffer 0.136100 -0.2827 0.02861
## 250 meter buffer 0.109200 -0.2471 -0.04877
## cnd(scl(trl)) cnd(scl(wll)) df logLik AICc delta weight
## 1500 meter buffer 0.16610 0.34340 9 -565.155 1149.1 0.00 0.732
## 1250 meter buffer 0.14970 0.32410 9 -566.478 1151.8 2.64 0.195
## 1750 meter buffer 0.17130 0.34550 9 -567.722 1154.3 5.13 0.056
## 4000 meter buffer 0.23450 0.42800 9 -570.437 1159.7 10.56 0.004
## 3750 meter buffer 0.23610 0.42510 9 -570.442 1159.7 10.57 0.004
## 2000 meter buffer 0.17170 0.35990 9 -570.687 1160.2 11.06 0.003
## 2250 meter buffer 0.19010 0.38950 9 -571.249 1161.3 12.19 0.002
## 4250 meter buffer 0.22410 0.42550 9 -571.525 1161.9 12.74 0.001
## 3500 meter buffer 0.23060 0.41520 9 -571.732 1162.3 13.15 0.001
## 3250 meter buffer 0.22160 0.41500 9 -572.507 1163.8 14.70 0.000
## 4500 meter buffer 0.21740 0.42380 9 -572.561 1163.9 14.81 0.000
## 2500 meter buffer 0.20200 0.39440 9 -572.805 1164.4 15.30 0.000
## 2750 meter buffer 0.20050 0.41580 9 -573.186 1165.2 16.06 0.000
## 1000 meter buffer 0.14490 0.29200 9 -573.203 1165.2 16.09 0.000
## 3000 meter buffer 0.20430 0.41690 9 -573.359 1165.5 16.41 0.000
## 4750 meter buffer 0.21630 0.42010 9 -574.034 1166.9 17.76 0.000
## 5000 meter buffer 0.21750 0.43480 9 -574.372 1167.6 18.43 0.000
## 750 meter buffer 0.08336 0.27470 9 -584.373 1187.6 38.44 0.000
## 500 meter buffer 0.09759 0.25800 9 -590.410 1199.6 50.51 0.000
## 250 meter buffer 0.18620 0.06643 9 -598.589 1216.0 66.87 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(deer_anthro_mods$`1500 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~ scale(harvest) +
## scale(seismic_lines) + scale(seismic_lines_3D) + scale(trails) +
## scale(wells) + scale(osm_industrial) + scale(pipeline_transmission_lines) +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 1148.3 1179.3 -565.2 1130.3 223
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 1.779 1.334
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.19497 0.54694 -0.356 0.72149
## scale(harvest) -0.08754 0.06511 -1.344 0.17879
## scale(seismic_lines) -0.07579 0.06319 -1.199 0.23039
## scale(seismic_lines_3D) -0.42869 0.09018 -4.754 2.00e-06 ***
## scale(trails) 0.16605 0.05173 3.210 0.00133 **
## scale(wells) 0.34337 0.07246 4.739 2.15e-06 ***
## scale(osm_industrial) 0.32669 0.05488 5.953 2.63e-09 ***
## scale(pipeline_transmission_lines) 0.06892 0.07549 0.913 0.36126
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Landscape
Now we will duplicate this analysis but for data that only includes landscape features.
For now we have to select EITHER lc_forest or lc_shrub as they are highly correlated
Black bear
bear_land_mods <- osm_landscape_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
# Landscape classes that aren't correlated
scale(lc_forest) +
scale(lc_grassland) +
scale(lc_developed) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(bear_land_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(lc_dvl)) cnd(scl(lc_frs))
## 250 meter buffer -0.5940 + -0.22460 -0.166500
## 500 meter buffer -0.5921 + -0.12600 -0.062020
## 1000 meter buffer -0.5903 + -0.08993 0.008256
## 1250 meter buffer -0.5899 + -0.09248 -0.009789
## 750 meter buffer -0.5913 + -0.10330 -0.005931
## 1500 meter buffer -0.5893 + -0.08210 -0.007234
## 2000 meter buffer -0.5886 + -0.06772 0.020710
## 2250 meter buffer -0.5881 + -0.06176 0.020950
## 2500 meter buffer -0.5883 + -0.06529 0.017800
## 1750 meter buffer -0.5889 + -0.06971 0.008760
## 2750 meter buffer -0.5889 + -0.07668 0.014030
## 3000 meter buffer -0.5890 + -0.07806 0.008553
## 3500 meter buffer -0.5892 + -0.08398 -0.015140
## 3750 meter buffer -0.5895 + -0.08584 -0.024250
## 3250 meter buffer -0.5890 + -0.07955 -0.004495
## 4000 meter buffer -0.5898 + -0.08630 -0.028750
## 4250 meter buffer -0.5900 + -0.08116 -0.034460
## 4500 meter buffer -0.5901 + -0.07535 -0.042560
## 5000 meter buffer -0.5897 + -0.06409 -0.055880
## 4750 meter buffer -0.5896 + -0.06502 -0.048360
## cnd(scl(lc_grs)) df logLik AICc delta weight
## 250 meter buffer -0.0338900 5 -456.395 923.1 0.00 0.915
## 500 meter buffer 0.0189700 5 -460.990 932.2 9.19 0.009
## 1000 meter buffer -0.0321200 5 -461.081 932.4 9.37 0.008
## 1250 meter buffer -0.0352000 5 -461.146 932.6 9.50 0.008
## 750 meter buffer -0.0007373 5 -461.415 933.1 10.04 0.006
## 1500 meter buffer -0.0284800 5 -461.545 933.4 10.30 0.005
## 2000 meter buffer -0.0290800 5 -461.624 933.5 10.46 0.005
## 2250 meter buffer -0.0354900 5 -461.691 933.6 10.59 0.005
## 2500 meter buffer -0.0274600 5 -461.799 933.9 10.81 0.004
## 1750 meter buffer -0.0223100 5 -461.833 933.9 10.88 0.004
## 2750 meter buffer -0.0054350 5 -461.882 934.0 10.97 0.004
## 3000 meter buffer -0.0036750 5 -461.942 934.1 11.09 0.004
## 3500 meter buffer 0.0004516 5 -462.001 934.3 11.21 0.003
## 3750 meter buffer 0.0057620 5 -462.011 934.3 11.23 0.003
## 3250 meter buffer -0.0014420 5 -462.046 934.4 11.30 0.003
## 4000 meter buffer 0.0147900 5 -462.049 934.4 11.31 0.003
## 4250 meter buffer 0.0192800 5 -462.178 934.6 11.57 0.003
## 4500 meter buffer 0.0195800 5 -462.252 934.8 11.71 0.003
## 5000 meter buffer 0.0020860 5 -462.255 934.8 11.72 0.003
## 4750 meter buffer 0.0046210 5 -462.334 934.9 11.88 0.002
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(bear_land_mods$`250 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(black_bear, absent_black_bear) ~ scale(lc_forest) + scale(lc_grassland) +
## scale(lc_developed) + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 922.8 940.0 -456.4 912.8 227
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.1092 0.3304
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.59401 0.14308 -4.151 3.3e-05 ***
## scale(lc_forest) -0.16655 0.06698 -2.487 0.012898 *
## scale(lc_grassland) -0.03389 0.05487 -0.618 0.536799
## scale(lc_developed) -0.22456 0.06436 -3.489 0.000485 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Caribou
caribou_land_mods <- osm_landscape_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(caribou, absent_caribou) ~
# Landscape classes that aren't correlated
scale(lc_forest) +
scale(lc_grassland) +
scale(lc_developed) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(caribou_land_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(lc_dvl)) cnd(scl(lc_frs))
## 250 meter buffer -5.040 + -0.07921 0.221800
## 1250 meter buffer -5.073 + -0.47570 0.118300
## 1500 meter buffer -5.073 + -0.46440 0.121600
## 500 meter buffer -5.021 + 0.02154 0.240000
## 750 meter buffer -5.023 + -0.10810 0.197600
## 1750 meter buffer -5.074 + -0.52210 0.106200
## 1000 meter buffer -5.039 + -0.33120 0.142100
## 2000 meter buffer -5.036 + -0.41610 0.090470
## 4750 meter buffer -5.016 + 0.02212 -0.050300
## 4500 meter buffer -5.022 + 0.09100 -0.038890
## 2250 meter buffer -5.017 + -0.33760 0.074890
## 5000 meter buffer -5.021 + 0.09191 -0.059830
## 4250 meter buffer -5.018 + 0.10870 -0.035480
## 2500 meter buffer -5.007 + -0.29410 0.057850
## 2750 meter buffer -5.002 + -0.26290 0.043150
## 4000 meter buffer -5.015 + 0.12760 -0.038380
## 3750 meter buffer -4.996 + -0.01350 -0.023500
## 3500 meter buffer -4.993 + -0.05918 -0.001797
## 3000 meter buffer -4.993 + -0.19810 0.032110
## 3250 meter buffer -4.988 + -0.10340 0.015120
## cnd(scl(lc_grs)) df logLik AICc delta weight
## 250 meter buffer -0.15410 5 -150.585 311.4 0.00 0.222
## 1250 meter buffer -0.12060 5 -151.199 312.7 1.23 0.120
## 1500 meter buffer -0.14280 5 -151.209 312.7 1.25 0.119
## 500 meter buffer -0.12290 5 -151.388 313.0 1.61 0.099
## 750 meter buffer -0.15220 5 -151.566 313.4 1.96 0.083
## 1750 meter buffer -0.10120 5 -151.644 313.6 2.12 0.077
## 1000 meter buffer -0.09792 5 -151.919 314.1 2.67 0.058
## 2000 meter buffer -0.08011 5 -152.577 315.4 3.98 0.030
## 4750 meter buffer -0.20900 5 -152.915 316.1 4.66 0.022
## 4500 meter buffer -0.21200 5 -152.952 316.2 4.73 0.021
## 2250 meter buffer -0.08965 5 -152.981 316.2 4.79 0.020
## 5000 meter buffer -0.19530 5 -153.011 316.3 4.85 0.020
## 4250 meter buffer -0.18810 5 -153.133 316.5 5.10 0.017
## 2500 meter buffer -0.09180 5 -153.206 316.7 5.24 0.016
## 2750 meter buffer -0.09703 5 -153.316 316.9 5.46 0.014
## 4000 meter buffer -0.15640 5 -153.335 316.9 5.50 0.014
## 3750 meter buffer -0.13740 5 -153.455 317.2 5.74 0.013
## 3500 meter buffer -0.13120 5 -153.480 317.2 5.79 0.012
## 3000 meter buffer -0.09639 5 -153.485 317.2 5.80 0.012
## 3250 meter buffer -0.10320 5 -153.602 317.5 6.03 0.011
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(caribou_land_mods$`250 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(caribou, absent_caribou) ~ scale(lc_forest) + scale(lc_grassland) +
## scale(lc_developed) + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 311.2 328.4 -150.6 301.2 227
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 5.351 2.313
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.03974 1.10127 -4.576 4.73e-06 ***
## scale(lc_forest) 0.22177 0.15873 1.397 0.162
## scale(lc_grassland) -0.15412 0.16528 -0.932 0.351
## scale(lc_developed) -0.07921 0.20030 -0.395 0.693
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Coyote
coyote_land_mods <- osm_landscape_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(coyote, absent_coyote) ~
# Landscape classes that aren't correlated
scale(lc_forest) +
scale(lc_grassland) +
scale(lc_developed) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(coyote_land_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(lc_dvl)) cnd(scl(lc_frs))
## 3750 meter buffer -1.369 + 0.4530 -0.2505
## 4500 meter buffer -1.366 + 0.4741 -0.2354
## 4750 meter buffer -1.365 + 0.4811 -0.2306
## 5000 meter buffer -1.364 + 0.4857 -0.2255
## 3500 meter buffer -1.368 + 0.4484 -0.2468
## 4250 meter buffer -1.367 + 0.4594 -0.2404
## 4000 meter buffer -1.369 + 0.4494 -0.2481
## 3250 meter buffer -1.368 + 0.4331 -0.2468
## 3000 meter buffer -1.368 + 0.4248 -0.2441
## 1250 meter buffer -1.367 + 0.4173 -0.1826
## 1750 meter buffer -1.364 + 0.4390 -0.1670
## 2750 meter buffer -1.366 + 0.4197 -0.2305
## 2000 meter buffer -1.363 + 0.4372 -0.1757
## 1500 meter buffer -1.365 + 0.4236 -0.1716
## 2250 meter buffer -1.364 + 0.4356 -0.1913
## 2500 meter buffer -1.365 + 0.4240 -0.2121
## 1000 meter buffer -1.369 + 0.3838 -0.1960
## 750 meter buffer -1.370 + 0.3340 -0.1862
## 500 meter buffer -1.370 + 0.2110 -0.1916
## 250 meter buffer -1.373 + 0.1435 -0.1977
## cnd(scl(lc_grs)) df logLik AICc delta weight
## 3750 meter buffer -0.075340 5 -486.114 982.5 0.00 0.192
## 4500 meter buffer -0.092240 5 -486.175 982.6 0.12 0.181
## 4750 meter buffer -0.090670 5 -486.181 982.6 0.14 0.180
## 5000 meter buffer -0.091480 5 -486.459 983.2 0.69 0.136
## 3500 meter buffer -0.068580 5 -486.708 983.7 1.19 0.106
## 4250 meter buffer -0.083470 5 -486.996 984.3 1.76 0.080
## 4000 meter buffer -0.074200 5 -487.011 984.3 1.80 0.078
## 3250 meter buffer -0.052270 5 -488.095 986.5 3.96 0.027
## 3000 meter buffer -0.050350 5 -489.189 988.6 6.15 0.009
## 1250 meter buffer -0.017410 5 -490.437 991.1 8.65 0.003
## 1750 meter buffer -0.050530 5 -490.523 991.3 8.82 0.002
## 2750 meter buffer -0.055640 5 -490.938 992.1 9.65 0.002
## 2000 meter buffer -0.072880 5 -491.295 992.9 10.36 0.001
## 1500 meter buffer -0.028570 5 -491.331 992.9 10.43 0.001
## 2250 meter buffer -0.083070 5 -491.414 993.1 10.60 0.001
## 2500 meter buffer -0.070990 5 -491.641 993.5 11.05 0.001
## 1000 meter buffer 0.009004 5 -491.801 993.9 11.38 0.001
## 750 meter buffer 0.048060 5 -497.968 1006.2 23.71 0.000
## 500 meter buffer 0.039140 5 -513.369 1037.0 54.51 0.000
## 250 meter buffer 0.006203 5 -519.650 1049.6 67.07 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(coyote_land_mods$`3750 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(coyote, absent_coyote) ~ scale(lc_forest) + scale(lc_grassland) +
## scale(lc_developed) + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 982.2 999.5 -486.1 972.2 227
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.1374 0.3707
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.36876 0.16060 -8.523 < 2e-16 ***
## scale(lc_forest) -0.25046 0.05751 -4.355 1.33e-05 ***
## scale(lc_grassland) -0.07534 0.06625 -1.137 0.255
## scale(lc_developed) 0.45296 0.06442 7.032 2.04e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fisher
fisher_land_mods <- osm_landscape_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(fisher, absent_fisher) ~
# Landscape classes that aren't correlated
scale(lc_forest) +
scale(lc_grassland) +
scale(lc_developed) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(fisher_land_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(lc_dvl)) cnd(scl(lc_frs))
## 2500 meter buffer -2.968 + 0.4462 0.4626
## 3250 meter buffer -2.965 + 0.4789 0.4209
## 2750 meter buffer -2.967 + 0.4574 0.4540
## 3000 meter buffer -2.965 + 0.4731 0.4321
## 2250 meter buffer -2.969 + 0.4208 0.4861
## 3500 meter buffer -2.960 + 0.4585 0.4082
## 2000 meter buffer -2.971 + 0.4007 0.5042
## 3750 meter buffer -2.952 + 0.4343 0.3842
## 4250 meter buffer -2.950 + 0.4292 0.3684
## 4000 meter buffer -2.950 + 0.4245 0.3701
## 4500 meter buffer -2.947 + 0.4115 0.3672
## 4750 meter buffer -2.946 + 0.3972 0.3681
## 1750 meter buffer -2.964 + 0.3642 0.4815
## 5000 meter buffer -2.944 + 0.3765 0.3687
## 1250 meter buffer -2.968 + 0.3575 0.5178
## 1000 meter buffer -2.981 + 0.3499 0.5883
## 750 meter buffer -2.988 + 0.3899 0.6380
## 500 meter buffer -2.988 + 0.4346 0.6501
## 1500 meter buffer -2.960 + 0.3511 0.4675
## 250 meter buffer -2.942 + 0.2841 0.4275
## cnd(scl(lc_grs)) df logLik AICc delta weight
## 2500 meter buffer -0.188800 5 -278.735 567.7 0.00 0.205
## 3250 meter buffer -0.264100 5 -278.875 568.0 0.28 0.178
## 2750 meter buffer -0.207000 5 -279.031 568.3 0.59 0.153
## 3000 meter buffer -0.244400 5 -279.062 568.4 0.66 0.148
## 2250 meter buffer -0.118400 5 -279.704 569.7 1.94 0.078
## 3500 meter buffer -0.248500 5 -279.715 569.7 1.96 0.077
## 2000 meter buffer -0.077940 5 -280.356 571.0 3.24 0.041
## 3750 meter buffer -0.220700 5 -280.825 571.9 4.18 0.025
## 4250 meter buffer -0.241300 5 -280.846 572.0 4.22 0.025
## 4000 meter buffer -0.226200 5 -281.074 572.4 4.68 0.020
## 4500 meter buffer -0.219800 5 -281.509 573.3 5.55 0.013
## 4750 meter buffer -0.198500 5 -281.828 573.9 6.19 0.009
## 1750 meter buffer -0.043210 5 -282.369 575.0 7.27 0.005
## 5000 meter buffer -0.167500 5 -282.423 575.1 7.38 0.005
## 1250 meter buffer 0.027420 5 -282.588 575.4 7.71 0.004
## 1000 meter buffer 0.079270 5 -282.699 575.7 7.93 0.004
## 750 meter buffer 0.159500 5 -282.743 575.8 8.02 0.004
## 500 meter buffer 0.248000 5 -282.965 576.2 8.46 0.003
## 1500 meter buffer -0.009224 5 -283.023 576.3 8.58 0.003
## 250 meter buffer 0.223600 5 -288.867 588.0 20.27 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(fisher_land_mods$`2500 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(fisher, absent_fisher) ~ scale(lc_forest) + scale(lc_grassland) +
## scale(lc_developed) + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 567.5 584.7 -278.7 557.5 227
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.5454 0.7385
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.96838 0.32100 -9.247 < 2e-16 ***
## scale(lc_forest) 0.46265 0.14544 3.181 0.00147 **
## scale(lc_grassland) -0.18883 0.12061 -1.566 0.11745
## scale(lc_developed) 0.44625 0.09245 4.827 1.38e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Grey wolf
wolf_land_mods <- osm_landscape_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(grey_wolf, absent_grey_wolf) ~
# Landscape classes that aren't correlated
scale(lc_forest) +
scale(lc_grassland) +
scale(lc_developed) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(wolf_land_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(lc_dvl)) cnd(scl(lc_frs))
## 3500 meter buffer -3.046 + -0.27920 0.005707
## 3750 meter buffer -3.046 + -0.30140 -0.014530
## 4000 meter buffer -3.042 + -0.33190 -0.025040
## 3250 meter buffer -3.034 + -0.26430 0.022660
## 4250 meter buffer -3.033 + -0.33770 -0.041090
## 4500 meter buffer -3.028 + -0.33980 -0.058000
## 3000 meter buffer -3.022 + -0.23720 0.021700
## 4750 meter buffer -3.020 + -0.30770 -0.064110
## 5000 meter buffer -3.015 + -0.28890 -0.079160
## 2750 meter buffer -3.012 + -0.20290 0.037410
## 2500 meter buffer -2.998 + -0.14740 0.055450
## 500 meter buffer -3.001 + -0.09243 -0.137000
## 2250 meter buffer -2.992 + -0.10050 0.088960
## 250 meter buffer -2.990 + -0.14610 -0.090580
## 2000 meter buffer -2.989 + -0.05724 0.105700
## 1750 meter buffer -2.986 + -0.01933 0.123700
## 1500 meter buffer -2.984 + 0.01408 0.116800
## 1250 meter buffer -2.982 + -0.02109 0.071710
## 750 meter buffer -2.985 + -0.07707 -0.018950
## 1000 meter buffer -2.983 + -0.03963 0.035650
## cnd(scl(lc_grs)) df logLik AICc delta weight
## 3500 meter buffer 0.353800 5 -259.379 529.0 0.00 0.243
## 3750 meter buffer 0.348300 5 -259.509 529.3 0.26 0.213
## 4000 meter buffer 0.332000 5 -259.717 529.7 0.68 0.173
## 3250 meter buffer 0.320400 5 -260.146 530.6 1.54 0.113
## 4250 meter buffer 0.298900 5 -260.386 531.0 2.01 0.089
## 4500 meter buffer 0.271000 5 -260.849 532.0 2.94 0.056
## 3000 meter buffer 0.276400 5 -261.257 532.8 3.76 0.037
## 4750 meter buffer 0.240300 5 -261.651 533.6 4.54 0.025
## 5000 meter buffer 0.216000 5 -262.105 534.5 5.45 0.016
## 2750 meter buffer 0.237400 5 -262.113 534.5 5.47 0.016
## 2500 meter buffer 0.172800 5 -263.343 537.0 7.93 0.005
## 500 meter buffer -0.167200 5 -263.612 537.5 8.47 0.004
## 2250 meter buffer 0.133200 5 -263.900 538.1 9.04 0.003
## 250 meter buffer -0.070320 5 -264.221 538.7 9.68 0.002
## 2000 meter buffer 0.108200 5 -264.240 538.7 9.72 0.002
## 1750 meter buffer 0.075180 5 -264.467 539.2 10.18 0.001
## 1500 meter buffer 0.029150 5 -264.645 539.6 10.53 0.001
## 1250 meter buffer 0.013030 5 -264.781 539.8 10.80 0.001
## 750 meter buffer 0.002586 5 -264.806 539.9 10.85 0.001
## 1000 meter buffer 0.003798 5 -264.843 540.0 10.93 0.001
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(wolf_land_mods$`3500 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(grey_wolf, absent_grey_wolf) ~ scale(lc_forest) + scale(lc_grassland) +
## scale(lc_developed) + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 528.8 546.0 -259.4 518.8 227
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.3526 0.5938
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.046161 0.267243 -11.398 < 2e-16 ***
## scale(lc_forest) 0.005707 0.149326 0.038 0.96951
## scale(lc_grassland) 0.353790 0.117829 3.003 0.00268 **
## scale(lc_developed) -0.279167 0.148037 -1.886 0.05932 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Lynx
lynx_land_mods <- osm_landscape_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(lynx, absent_lynx) ~
# Landscape classes that aren't correlated
scale(lc_forest) +
scale(lc_grassland) +
scale(lc_developed) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(lynx_land_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(lc_dvl)) cnd(scl(lc_frs))
## 500 meter buffer -1.859 + -0.10280 -0.0125600
## 750 meter buffer -1.858 + -0.06674 0.0417200
## 250 meter buffer -1.856 + -0.12880 -0.0552600
## 1000 meter buffer -1.856 + -0.04118 0.0430100
## 1250 meter buffer -1.854 + 0.01610 0.0536600
## 1500 meter buffer -1.854 + 0.04787 0.0516600
## 4750 meter buffer -1.851 + 0.04170 -0.0215400
## 5000 meter buffer -1.851 + 0.03511 -0.0158500
## 3000 meter buffer -1.855 + 0.01806 -0.0246500
## 4500 meter buffer -1.851 + 0.03534 -0.0272400
## 3250 meter buffer -1.855 + 0.01481 -0.0321000
## 1750 meter buffer -1.852 + 0.05850 0.0488100
## 2750 meter buffer -1.854 + 0.01892 -0.0171500
## 2000 meter buffer -1.852 + 0.05903 0.0424700
## 3500 meter buffer -1.854 + 0.02250 -0.0371500
## 3750 meter buffer -1.853 + 0.03749 -0.0357500
## 4000 meter buffer -1.852 + 0.03489 -0.0346700
## 4250 meter buffer -1.852 + 0.02934 -0.0315300
## 2250 meter buffer -1.852 + 0.05816 0.0232700
## 2500 meter buffer -1.853 + 0.03873 -0.0002094
## cnd(scl(lc_grs)) df logLik AICc delta weight
## 500 meter buffer 0.095240 5 -452.401 915.1 0.00 0.207
## 750 meter buffer 0.088070 5 -453.101 916.5 1.40 0.103
## 250 meter buffer 0.005896 5 -453.318 916.9 1.83 0.083
## 1000 meter buffer 0.078390 5 -453.571 917.4 2.34 0.064
## 1250 meter buffer 0.070890 5 -453.786 917.8 2.77 0.052
## 1500 meter buffer 0.051570 5 -453.932 918.1 3.06 0.045
## 4750 meter buffer -0.090990 5 -453.940 918.1 3.08 0.044
## 5000 meter buffer -0.092860 5 -453.973 918.2 3.15 0.043
## 3000 meter buffer 0.050080 5 -454.185 918.6 3.57 0.035
## 4500 meter buffer -0.066000 5 -454.204 918.7 3.61 0.034
## 3250 meter buffer 0.045350 5 -454.221 918.7 3.64 0.033
## 1750 meter buffer 0.011150 5 -454.270 918.8 3.74 0.032
## 2750 meter buffer 0.044980 5 -454.295 918.9 3.79 0.031
## 2000 meter buffer 0.002748 5 -454.338 918.9 3.88 0.030
## 3500 meter buffer 0.023660 5 -454.340 918.9 3.88 0.030
## 3750 meter buffer -0.006762 5 -454.377 919.0 3.95 0.029
## 4000 meter buffer -0.025900 5 -454.396 919.1 3.99 0.028
## 4250 meter buffer -0.035690 5 -454.417 919.1 4.03 0.028
## 2250 meter buffer -0.006766 5 -454.448 919.2 4.09 0.027
## 2500 meter buffer 0.014420 5 -454.526 919.3 4.25 0.025
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(lynx_land_mods$`500 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(lynx, absent_lynx) ~ scale(lc_forest) + scale(lc_grassland) +
## scale(lc_developed) + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 914.8 932.0 -452.4 904.8 227
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.1629 0.4037
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.85919 0.17480 -10.636 <2e-16 ***
## scale(lc_forest) -0.01256 0.06974 -0.180 0.857
## scale(lc_grassland) 0.09524 0.05849 1.628 0.103
## scale(lc_developed) -0.10278 0.08084 -1.271 0.204
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Moose
moose_land_mods <- osm_landscape_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(moose, absent_moose) ~
# Landscape classes that aren't correlated
scale(lc_forest) +
scale(lc_grassland) +
scale(lc_developed) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(moose_land_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(lc_dvl)) cnd(scl(lc_frs))
## 250 meter buffer -1.824 + -0.38910 -0.1894000
## 500 meter buffer -1.811 + -0.26520 -0.1252000
## 1250 meter buffer -1.812 + -0.24850 -0.1302000
## 1000 meter buffer -1.809 + -0.24500 -0.0973000
## 750 meter buffer -1.807 + -0.22420 -0.0569800
## 1500 meter buffer -1.808 + -0.20490 -0.1309000
## 1750 meter buffer -1.808 + -0.20460 -0.1226000
## 2000 meter buffer -1.803 + -0.17250 -0.0978700
## 2250 meter buffer -1.803 + -0.17450 -0.0896700
## 5000 meter buffer -1.795 + -0.03136 0.0229600
## 2500 meter buffer -1.802 + -0.16470 -0.0705300
## 2750 meter buffer -1.801 + -0.15120 -0.0541200
## 4750 meter buffer -1.794 + -0.04654 0.0170200
## 3000 meter buffer -1.799 + -0.13010 -0.0322300
## 4500 meter buffer -1.794 + -0.06240 0.0152600
## 4250 meter buffer -1.794 + -0.06813 0.0144000
## 4000 meter buffer -1.795 + -0.07760 0.0091480
## 3250 meter buffer -1.797 + -0.11400 -0.0159600
## 3500 meter buffer -1.796 + -0.10190 -0.0070070
## 3750 meter buffer -1.795 + -0.08647 -0.0006682
## cnd(scl(lc_grs)) df logLik AICc delta weight
## 250 meter buffer -0.084560 5 -435.539 881.3 0.00 0.988
## 500 meter buffer 0.012840 5 -441.178 892.6 11.28 0.004
## 1250 meter buffer 0.045990 5 -441.280 892.8 11.48 0.003
## 1000 meter buffer 0.019960 5 -441.850 894.0 12.62 0.002
## 750 meter buffer 0.032590 5 -442.554 895.4 14.03 0.001
## 1500 meter buffer 0.041120 5 -442.573 895.4 14.07 0.001
## 1750 meter buffer 0.048510 5 -442.699 895.7 14.32 0.001
## 2000 meter buffer 0.023710 5 -443.902 898.1 16.73 0.000
## 2250 meter buffer 0.013000 5 -443.941 898.1 16.80 0.000
## 5000 meter buffer -0.150100 5 -444.265 898.8 17.45 0.000
## 2500 meter buffer 0.020290 5 -444.293 898.9 17.51 0.000
## 2750 meter buffer 0.037060 5 -444.686 899.6 18.29 0.000
## 4750 meter buffer -0.108600 5 -444.907 900.1 18.74 0.000
## 3000 meter buffer 0.022790 5 -445.144 900.6 19.21 0.000
## 4500 meter buffer -0.078030 5 -445.167 900.6 19.26 0.000
## 4250 meter buffer -0.067810 5 -445.207 900.7 19.34 0.000
## 4000 meter buffer -0.049860 5 -445.328 900.9 19.58 0.000
## 3250 meter buffer 0.012890 5 -445.383 901.0 19.69 0.000
## 3500 meter buffer -0.002757 5 -445.478 901.2 19.88 0.000
## 3750 meter buffer -0.022080 5 -445.552 901.4 20.03 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(moose_land_mods$`250 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(moose, absent_moose) ~ scale(lc_forest) + scale(lc_grassland) +
## scale(lc_developed) + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 881.1 898.3 -435.5 871.1 227
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.2284 0.4779
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.82437 0.20458 -8.917 < 2e-16 ***
## scale(lc_forest) -0.18942 0.07609 -2.489 0.0128 *
## scale(lc_grassland) -0.08456 0.06399 -1.321 0.1864
## scale(lc_developed) -0.38907 0.08615 -4.516 6.3e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Red fox
fox_land_mods <- osm_landscape_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(red_fox, absent_red_fox) ~
# Landscape classes that aren't correlated
scale(lc_forest) +
scale(lc_grassland) +
scale(lc_developed) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(fox_land_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(lc_dvl)) cnd(scl(lc_frs))
## 4750 meter buffer -4.312 + 0.2866 -0.47070
## 5000 meter buffer -4.314 + 0.2946 -0.49010
## 4500 meter buffer -4.306 + 0.2758 -0.46010
## 4250 meter buffer -4.289 + 0.2666 -0.43620
## 3750 meter buffer -4.268 + 0.2492 -0.31340
## 4000 meter buffer -4.265 + 0.2689 -0.37450
## 3500 meter buffer -4.258 + 0.2469 -0.28260
## 3250 meter buffer -4.259 + 0.2281 -0.24610
## 1500 meter buffer -4.286 + 0.1721 -0.14490
## 2750 meter buffer -4.251 + 0.2309 -0.25520
## 2500 meter buffer -4.280 + 0.2042 -0.21590
## 2250 meter buffer -4.293 + 0.2030 -0.17680
## 3000 meter buffer -4.239 + 0.2446 -0.24900
## 1750 meter buffer -4.291 + 0.1625 -0.18910
## 1250 meter buffer -4.243 + 0.2205 -0.14630
## 2000 meter buffer -4.283 + 0.2067 -0.17390
## 1000 meter buffer -4.226 + 0.2698 -0.11050
## 750 meter buffer -4.198 + 0.3321 0.04518
## 500 meter buffer -4.126 + 0.6055 0.25200
## 250 meter buffer -4.126 + 0.4755 0.02398
## cnd(scl(lc_grs)) df logLik AICc delta weight
## 4750 meter buffer 0.4441 5 -172.612 355.5 0.00 0.325
## 5000 meter buffer 0.4375 5 -172.740 355.7 0.26 0.286
## 4500 meter buffer 0.4370 5 -173.157 356.6 1.09 0.188
## 4250 meter buffer 0.4273 5 -174.117 358.5 3.01 0.072
## 3750 meter buffer 0.4867 5 -174.819 359.9 4.41 0.036
## 4000 meter buffer 0.4291 5 -175.095 360.5 4.97 0.027
## 3500 meter buffer 0.4829 5 -175.667 361.6 6.11 0.015
## 3250 meter buffer 0.5048 5 -176.028 362.3 6.83 0.011
## 1500 meter buffer 0.5401 5 -176.163 362.6 7.10 0.009
## 2750 meter buffer 0.4595 5 -176.645 363.6 8.07 0.006
## 2500 meter buffer 0.5171 5 -176.796 363.9 8.37 0.005
## 2250 meter buffer 0.5380 5 -176.816 363.9 8.41 0.005
## 3000 meter buffer 0.4499 5 -176.831 363.9 8.44 0.005
## 1750 meter buffer 0.5305 5 -177.248 364.8 9.27 0.003
## 1250 meter buffer 0.4693 5 -177.375 365.0 9.53 0.003
## 2000 meter buffer 0.5172 5 -177.500 365.3 9.78 0.002
## 1000 meter buffer 0.4439 5 -178.281 366.8 11.34 0.001
## 750 meter buffer 0.4519 5 -178.480 367.2 11.74 0.001
## 500 meter buffer 0.3413 5 -180.448 371.2 15.67 0.000
## 250 meter buffer 0.2554 5 -181.329 372.9 17.43 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)summary(fox_land_mods$`4750 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(red_fox, absent_red_fox) ~ scale(lc_forest) + scale(lc_grassland) +
## scale(lc_developed) + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 355.2 372.5 -172.6 345.2 227
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 1.299 1.14
## Number of obs: 232, groups: array, 6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.3117 0.5237 -8.233 <2e-16 ***
## scale(lc_forest) -0.4707 0.1928 -2.441 0.0147 *
## scale(lc_grassland) 0.4441 0.1794 2.476 0.0133 *
## scale(lc_developed) 0.2866 0.1200 2.388 0.0170 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1White-tailed deer
deer_land_mods <- osm_landscape_df_2021_2022 %>%
# use purrr map to run the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~
# Landscape classes that aren't correlated
scale(lc_forest) +
scale(lc_grassland) +
scale(lc_developed) +
# Random effect of array
(1|array),
data = .,
family = 'binomial'))
# run model selection
model.sel(deer_land_mods)## Model selection table
## cnd((Int)) dsp((Int)) cnd(scl(lc_dvl)) cnd(scl(lc_frs))
## 4500 meter buffer -0.1624 + 0.20190 -0.20540
## 4250 meter buffer -0.1625 + 0.19590 -0.19220
## 4750 meter buffer -0.1620 + 0.22300 -0.21080
## 4000 meter buffer -0.1635 + 0.19380 -0.17740
## 5000 meter buffer -0.1623 + 0.23070 -0.22070
## 3750 meter buffer -0.1641 + 0.19650 -0.16600
## 3500 meter buffer -0.1642 + 0.19470 -0.15780
## 3250 meter buffer -0.1642 + 0.19220 -0.15230
## 3000 meter buffer -0.1632 + 0.19470 -0.14750
## 2750 meter buffer -0.1626 + 0.19930 -0.14240
## 2500 meter buffer -0.1626 + 0.19860 -0.13980
## 2250 meter buffer -0.1623 + 0.20390 -0.13610
## 2000 meter buffer -0.1618 + 0.19850 -0.12930
## 1750 meter buffer -0.1616 + 0.18840 -0.12270
## 1500 meter buffer -0.1628 + 0.16270 -0.12320
## 1250 meter buffer -0.1616 + 0.16220 -0.10650
## 1000 meter buffer -0.1599 + 0.16350 -0.07477
## 750 meter buffer -0.1615 + 0.10970 -0.08782
## 250 meter buffer -0.1670 + -0.04939 -0.18920
## 500 meter buffer -0.1652 + 0.01760 -0.13120
## cnd(scl(lc_grs)) df logLik AICc delta weight
## 4500 meter buffer 0.37070 5 -587.931 1186.1 0.00 0.611
## 4250 meter buffer 0.36210 5 -589.295 1188.9 2.73 0.156
## 4750 meter buffer 0.34440 5 -589.330 1188.9 2.80 0.151
## 4000 meter buffer 0.35440 5 -590.591 1191.4 5.32 0.043
## 5000 meter buffer 0.31720 5 -591.140 1192.5 6.42 0.025
## 3750 meter buffer 0.34090 5 -592.009 1194.3 8.16 0.010
## 3500 meter buffer 0.32930 5 -593.436 1197.1 11.01 0.002
## 3250 meter buffer 0.32510 5 -593.863 1198.0 11.86 0.002
## 3000 meter buffer 0.30920 5 -595.470 1201.2 15.08 0.000
## 2750 meter buffer 0.28810 5 -597.459 1205.2 19.06 0.000
## 2500 meter buffer 0.27620 5 -599.026 1208.3 22.19 0.000
## 2250 meter buffer 0.25030 5 -601.146 1212.6 26.43 0.000
## 2000 meter buffer 0.23720 5 -602.343 1215.0 28.82 0.000
## 1750 meter buffer 0.22310 5 -603.842 1218.0 31.82 0.000
## 1500 meter buffer 0.21270 5 -606.062 1222.4 36.26 0.000
## 1250 meter buffer 0.19350 5 -608.795 1227.9 41.73 0.000
## 1000 meter buffer 0.14700 5 -614.226 1238.7 52.59 0.000
## 750 meter buffer 0.12650 5 -617.400 1245.1 58.94 0.000
## 250 meter buffer 0.03774 5 -618.659 1247.6 61.46 0.000
## 500 meter buffer 0.08188 5 -620.265 1250.8 64.67 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)In summary it does matter whether you are looking at landscape or anthropogenic features which spatial scale you use. Only two species (lynx and moose) had spatial scales for both top models within 1000m (1km). There were also more well-defined top models (delta AIC greater than 2.00) when you subset data to specifically look at anthro or landscape features.